Defense  Nuclear  Agency 
Alexandria,  VA  22310-3398 


DNA-TR-93-187 


Material  Modeling  in  the  CRALE  Code 

Users’  Manual 


Shel  Schuster 

The  Titan  Corporation 

Titan  Research  &  Technology  Division 

9410  Topanga  Canyon  Blvd,  Suite  104 

Chatsworth,  CA  91311-5771 


December  1994 


Technical  Report 


CONTRACT  No.  DNA  001 -93-C-0055 

f 


Approved  for  public  release; 
distribution  is  unlimited. 


19941219  071 


Destroy  this  report  when  it  is  no  longer  needed.  Do  not 
return  to  sender. 


PLEASE  NOTIFY  THE  DEFENSE  NUCLEAR  AGENCY. 
ATTN:  CSTI,  6801  TELEGRAPH  ROAD.  ALEXANDRIA.  VA 
22310-3398.  IF  YOUR  ADDRESS  IS  INCORRECT.  IF  YOU 
WISH  IT  DELETED  FROM  THE  DISTRIBUTION  LIST.  OR 
IFTHE  ADDRESSEE  IS  NO  LONGER  EMPLOYED  BY  YOUR 
ORGANIZATION. 


CUT  HERE  AND  RETURN 


DISTRIBUTION  LIST  UPDATE 


This  mailer  is  provided  to  enable  DNA  to  maintain  current  distribution  lists  for  reports.  (We  would 
appreciate  your  providing  the  requested  information.) 


□  Add  the  individual  listed  to  your  distribution  list. 

□  Delete  the  cited  organization /individual. 

□  Change  of  address. 


NAME:  - - - - - — - - 

ORGANIZATION:  - - - - - - - 

OLD  ADDRESS  CURRENT  ADDRESS 


'  TELEPHONE  NUMBER:  1 - 1 


DNA  PUBLICATION  NUMBER/TITLE  C H AN G ES/DELETION S/ AD D ITIO NS,  etc.) 


DNA  OR  OTHER  GOVERNMENT  CONTRACT  NUMBER:  - - - 

CERTIFICATION  OF  NEED-TO-KNOW  BY  GOVERNMENT  SPONSOR  (if  other  than  DNA): 
|  SPONSORING  ORGANIZATION:  - - - - - - 

I 

CONTRACTING  OFFICER  OR  REPRESENTATIVE:  - - 


NOTE: 

Please  return  the  mailing  label  from  the 
document  so  that  any  additions,  changes, 
corrections  or  deletions  can  be  made  easily. 
For  distribution  cancellation  or  more 
information  call  DNA/IMAS  (703)  325-1036. 


SIGNATURE: 


DEFENSE  NUCLEAR  AGENCY 
ATTN:  I  MAS 

6801  TELEGRAPH  ROAD 
ALEXANDRIA,  VA  22310-3398 


DEFENSE  NUCLEAR  AGENCY 
ATTN:  IMAS 

6801  TELEGRAPH  ROAD 
ALEXANDRIA,  VA  22310-3398 


REPORT  DOCUMENTATION  PAGE  Townee 

'  Pubkc  .WnK.s  bu.»n  to  th.  coll.B.on  d  format, or  .  1C  .v.rag.  1  hour  per  -espons.  .«clucJ.«8  '»•  !«"•  'o'  -.v»..rc  ■'’*<'»««>"*  »“'<*'"« 

j  th*  data  needed  and  omptetine  and  f*v*w.ng  th*  colk*d.on  ol  information  S*nd  comments  regarding  this  buroen  est.mat*  ot  any  other  asp*d  o  tfus  co 

StZSTSJS. *  -ZSe  .H»  burden,  1oW«h,n9lon  Headquarters  Screes  D...c,o,ar.  to  tn«>.<h.„on  Ope.ato, ■  «- *  *£"*  ««  D™  H^' 

1204.  Arlington.  VA  22202-4302.  and  1o  the  Otfce  of  Management  and  Budget,  Paperwork  Reduction  Project  {0704-0166),  Washington  DC  20503  _ _ _ 

1.  AGENCY  USE  ONLY  [Leave  bilrt)  2.  REPORT  DATE  1  REPORT  TYPE  AND  DATES  COVERED 

I  941201  Technical  930312  -  931212 


4.  TITLE  AND  SUBTITLE 

Material  Modeling  in  the  CRALE  Code 
Users'  Manual 


6.  AUTHOR(S) 

Shel  Schuster 


5.  FUNDING  NUMBERS 

c  -  DNA  001-93-C-0055 
PE  -  62715H 


PR  -  CD,  CD,  CD 
TA  -  CD,  CC,  CD 
WU  -  DH600200 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

The  Titan_  Corporation 

Titan  Research  &  Technology  Division 
9410  Topanga  Canyon  Blvd,  Suite  104 
Chatsworth,  CA  91311-5771 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Defense  Nuclear  Agency 
6801  Telegraph  Road 
Alexandria,  VA  22310-3398 
FCTT/Rinehart 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

TRT  3336TR-1 


*0.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


DNA-TR-93-187 


n.  supplementary  NOTES  This  wQrk  wag  sponsored  by  the  Defense  Nuclear  Agency 

under  RDT&E  RMC  Codes  T4613D  CD  CD  60020  5900A  25904D,  T4662D  CD 
CC  60405  572 0A  25904D  and  T4613D  CD  CD  60027  5900A  25904D. 

12a.  DISTRIBUTION/AVAILABILITY  STATEMENT  "  "  l2b-  DISTRIBUTION  CODE 

Approved  for  public  release;  distribution  is 
unlimited . 


13.  ABSTFIACT  [Maximum  200  words) 

Finite  difference  for  continuum  mechanics  codes  are  capable  of 
providing  solutions  for  the  behavior  of  a  wide  range  of  materials  under 
a  variety  of  boundary  and  initial  conditions.  The  key  to  obtaining 
accurate  solutions  lies  in  the  fidelity  of  models  used  to  describe  the 
material's  response  to  the  strains  imposed,  i.e.,  its  Equation  of  State 
(EOS) .  This  report  presents  both  the  general  EOS  alqorithms  used  m  the 
CRALE  1-  and  2D  codes  and  several  specific  EOS  models  currently  in  use. 

A  driver  routine  that  exercises  these  models  over  specific  load/unload 
paths  is  available  to  assist  in  developing  and  checking  the  parameters 
needed  for  materials  of  interest.  A  discussion  of  this  driver  and  a 
users'  manual  to  assist  in  its  use  are  also  presented. 

In  the  course  of  developing  EOS  models  of  interest  to  DNA  and  as 
part  of  the  TRT  participation  in  the  DNA  HYDROPLUS  program,  a  number  of 
data  bases  have  been  developed.  These  include  files  containing; 


IS.  NUMBER  OF  PAGES 
102 


16.  PRICE  CODE 


20.  LIMITATION  OF 
ABSTRACT 


Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Sta  239-16 
256-102 


14.  SUBJECT  TERMS 

CRALE 

EOS 

MATERIAL  MODELING 

,  HYDROPLUS 

USER'S  MANUAL 

i 

17.'  SECURITY  CLASSIFICATION 

OF  REPORT 

18.  SECURITY  CLASSIFICATION 

OF  THIS  PAGE 

19.  SECURITY  CLASSIFICATION 

OF  ABSTRACT 

UNCLASSIFIED 

UNCLASSIFIED 

UNCLASSIFIED 

NSN  7540-01-280-5500 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


CLASSIFIED  BY; 

N/A 

since 

Unclassified. 

DECLASSIFY  ON: 

N/A 

since 

Unclassified. 

experimental  shock  Hugoniots  and  release  paths,  parameters  for  EOS 
models  used  in  the  HYDROPLUS  program,  and  the  nuclear  shot  experimental 
and  calculational  data  bases  developed  for  HYDROPLUS.  A  description 
of  the  various  data  bases  with  instructions  for  "their  access  is  in¬ 
cluded  in  this  report. 


Accesicn  For 

NTIS  CRAll 
!  DUG  TA3 


A-/  i 


i 


II 


SECURITY  CLASSIFICATION  OP  THIS  PAGE 

UNCLASSIFIED 


CONVERSION  TABLE 


Conversion  factors  for  U.S.  Customary  to  metric  (SI)  units  of  measurement 


MULTIPLY - »  BY - >  TO  GET 

TO  GET  < - BY  < - DIVIDE 


angstrom 

1.000  000XE-10 

meters  (m) 

atmosphere  (normal) 

1.013  25  X  E+2 

kilo  pascal  (kPa) 

bar 

1.000  000  X  E+2 

kilo  pascal  (kPa) 

barn 

1.000  000  X  E  -28 

meter2  (m2) 

british  thermal  unit  (thermochemical) 

1.054  350  XE +3 

joule  (J) 

calorie  (thermochemical) 

4.184  000 

joule  (J) 

cal  (thermochemical)/cm2 

4.184  000  XE -2 

mega  joule/m2(MJ/m2) 

curie 

3.700  000  X  E  +1 

*giga  becquerel  (GBq) 

degree  (angle) 

1.745  329  XE -2 

radian  (rad) 

degree  Fahrenheit 

tk  =  t°f +  459.67)/1.8 

degree  kelvin  (k) 

electron  volt 

1.602  19  X  E  -19 

joule  (J) 

erg 

1.000  000  XE -7 

joule  (J) 

erg/second 

1.000  000  XE -7 

watt  (W) 

foot 

3.0048  000  X  E  -1 

meter  (m) 

fout-pound-force 

1.355  818 

joule  (J) 

gallon  (U.S.  liquid) 

3.785  412  XE -3 

meter3  (m3) 

inch 

2.540  000  X  E  -2 

meter  (m) 

jerk 

1.000  000  XE +9 

joule  (J) 

joule/kilogram  (J/kg) 

1.000  000 

Gray  (gy) 

(radiation  dose  absorbed) 

kilotons 

4.183 

terajoules 

kip  (1000  Ibf) 

4.448  222  X  E  +3 

newton 

kip/inch2  (ksi) 

6.894  757  X  E  +3 

kilo  pascal  (kPa) 

ktap 

1.000  000  X  E+2 

newton-second/m2  (N-s/m2) 

micron 

1.000  000  XE -6 

meter  (m) 

mil 

2.540  000  X  E  -5 

meter  (m) 

mile  (international) 

1.609  344  XE +3 

meter  (m) 

ounce 

2.834  952  X  E  -2 

kilogram  (kg) 

pound-force  (lbs  avoirdupois) 

4.448  222 

newton  (N) 

pound-force  inch 

1.129  848XE-1 

newton-meter  (N-m) 

pound-force/inch 

1 .751  268  X  E  +2 

newton-meter  (N/m) 

pound-force/foot2 

4.788  026  X  E  -2 

kilo  pascal  (kPa) 

pound-force/inch2  (psi) 

6.894  747 

kilo  pascal  (kPa) 

pound-mass  (Ibm  avoirdupois) 

4.535  924  X  E  -1 

kilogram  (kg) 

pound-mass-foot2  (moment  of  inertia) 

4.214  011  XE-2 

kilogram-meter2  (kg-m2) 

pound-mass/foot3 

1.601  846XE+1 

kilogram/meter3  (kg-m3) 

rad  (radiation  dose  absorbed 

1.000  000  xE -2 

“Gray  (Gy) 

roentgen 

2.579  760  X  E  -4 

coulomb/kilogram  (* **C/Kg) 

shake 

1.000  000  XE -8 

seconds  (s) 

slug 

1.459  390XE+1 

kilogram  (kg) 

torr  (mm  Hg,  0°C) 

1.333  22XE-1 

kilo  pascal  (kPa) 

I 


I 


*  The  becquerel  (Bq)  is  the  SI  of  unit  of  radioactivity;  1  Bq  =  1  event/s. 

**  The  gray  (GY)  is  the  SI  unit  of  absorbed  radition. 


TABLE  OF  CONTENTS 


Section  Page 

Conversion  Table .  iii 

Figures .  vi 

Tables .  vii 

1  INTRODUCTION .  1 

1.1  Background .  1 

1 .2  Report  Structure .  2 

1 .3  Symbols  And  Definitions .  3 

2  GENERAL  EQUATION  OF  STATE  ALGORITHMS .  5 

2.1  Pressure .  5 

2.2  Elastic  Deviatoric  Behavior .  6 

2.3  Plastic  Yielding .  6 

2.4  Fracture .  7 

3  EOS  MODELS .  9 

3.1  Subroutine  Eqst:  Overview  And  Generic  Input .  9 

3.2  Specific  Eos  Models .  13 

3.2.1  TheSHELEOS .  15 

3.2.2  VIC .  21 

3.2.3  CIST .  22 

3.2.4  HE . . .  24 

3.2.5  MCST . 28 

3.2.6  SCUB .  28 

3.2.7  H20 .  28 

3.2.8  UEOS .  28 

3.2.9  CAP .  29 

3.2.10  MIX .  30 

3.2.11  Others .  30 


IV 


TABLE  OF  CONTENTS  (Continued) 


Section 


Page 


4 


5 


\ 


EQS  -  THE  EQUATION  OF  STATE  DRIVER 


4.1  Overview . 

4.2  Code  Modules 


31 

32 


4.2.1  EQS,  the  Main  Program . 

4.2.2  EQN . 

4.2.3  EXPHUG . 

4.2.4  FLIP  -  Reverse  Hugoniots . 

4.2.5  HUG  -  The  Hugoniot  Driver . 

4.2.6  REL  -  Release  Adiabats  from  the  Hugoniot. 

4.2.7  PATH  -  Special  Paths . 

4.2.8  PLOTN-Plotting  Options . 


32 

37 

39 

42 

44 

48 

52 

53 


DATA  BASES 


59 


5.1  Hugoniots . 

5.1.1  Experimental  Data  (ROCKHUG  and  METHUG). 

5.1.2  Shock  Velocity,  Us,  vs  Particle  Velocity, Up,  Fits 

(USUP.dat,  1-1-94) . 


5.2  Release  Adiabats  (Relad.dat,  1-1-94) 

5.3  Material  Models  (Eoss.dat)  1-1-94 . 

5.4  Hydroplus  Archival  Databases . 


5.4.1  Introduction . 

5.4.2  Equation  of  State/Numerical  Results 

5.4.3  Nuclear  Shot  Data . 


69 

69 

71 


Appendix 

A  LIST  OF  FILES  IN  UNCLASSIFIED  DIRECTORY . 

B  EXAMPLE  OF  THE  HYDROPLUS  NUCLEAR  SHOT  DATA  BASE 


I 

I 


V 


FIGURES 


Figure  Page 

3- 1  Release  Path  From  the  CJ  Point  of  TNT  Showing  the 

Contributions  of  the  Three  in  the  JWL  EOS .  26 

4- 1  Flowchart  of  EQS,  the  Equation  of  State  Driver  Routine .  33 

4-2  Example  of  Plot  Generated  by  ROUT=PLOT 

with  NP[1  ]=2  or  3 .  56 

4-3  Example  of  Plot  Generated  by  ROUT=PLOT  with  NP[1  ] — 5  ....  57 

4-4  Example  of  Plot  Generated  by  ROUT=PLOT  with  NP[1  ]— 6  ...  58 


I 

i 


VI 


TABLES 


Table  Pa9e 

1-1  Glossary  of  general  EOS  terms  and  symbols .  3 

3-1  Summary  of  input  for  the  EOS  models .  *1 9 

3-2  EOS  input  variables .  1 1 

3-3  General  parameters  used  by  all  equations-of-state .  12 

3-4  Summary  of  CK  block  parameters  of  current  EOS  models .  14 

3-5  Equation-of-state  parameters  for  SHEL  (EOS  1 ) .  20 

3-6  Equation-of-state  parameters  for  VIC  (EOS  2) .  21 

3-7  Equation-of-state  parameters  for  CIST  (EOS  3) .  23 

3-8  JWL  parameters  for  some  common  high  explosives .  27 

3-9  Equation-of-state  parameters  for  HE  (EOS  4) .  27 

3-10  Equation-of-state  parameters  for  CAP  (EOS  9) .  29 

3- 1 1  Equation-of-state  parameters  for  MIX  (EOS  10) .  30 

4- 1  Summary  of  the  input  to  EQS . 34 

4-2  Parameters  oh  line  4  of  EQS  input .  36 

4-3  Example  of  the  output  generated  by  ROUT  =  EQN .  41 

4- 4  Example  of  FLIP  printout  for  Borehole  Correction  Factors .  43 

5- 1  Shock  relationships .  60 

5-2  Materials  stored  in  the  ROCKHUG.dat  file  (Unit  10)  11-11-93 .  62 

5-3  Materials  stored  in  the  METHUG.dat  file  (Unit  11)  11-1 1-93 .  63 

5-4  Example  of  experimental  Hugoniot  data  stored 

on  the  ROCKHUG.dat  and  METHUG.dat  files .  64 

5-5  Material  Us-Up  Relations  Available  on  USUP.dat  File .  65 

i 

vii 


TABLES  (Continued) 


Table  Page 

5-6  Example  of  Data  on  RELAD.dat  (1-1-94) .  66 

5-7  List  of  105  Sets  of  Release  Adiabats  on  RELAD.dat  (1-1-94) .  67 

5-8  Listing  of  the  EOS  Models  in  the  EOSS.dat  Data  File, 

Tape  14  as  of  10/27/93 .  68 

5-9  Materials  in  HYDROPLUS  EOS  Database .  70 


VIII 


SECTION  1 


INTRODUCTION 


1.1  BACKGROUND. 

For  any  finite  difference  continuum  mechanics  code,  it  is  necessary  on  every 
cycle  to  compute  the  complete  state  of  stress  in  each  cell  in  order  to  generate  the 
forces  needed  to  advance  the  solution.  The  set  of  algorithms  used  to  calculate  these 
stresses  may  range  from  the  simplest  perfect-gas  law,  in  which  the  mean  stress  (i.e., 
pressure,  P)  is  expressed  as  a  function  of  only  two  independent  variables,  the  density 
(p)  and  specific  energy  (E),  to  one  that  computes  the  complete  nine-component  stress 
tensor  (ay)  as  a  function  of  the  energy,  the  current  strain  tensor  (ey),  and  the  path  taken 
from  the  material's  initial  state  to  its  present  one.  More  complex  models  may  even 
combine  several  material  components  and/or  phases  (e.g.,  soil  and  air  or  water),  use 
non-isotropic  stress-strain  relations,  or  incorporate  time-dependent  rate  effects,  as 
needed  by  the  specific  problem  at  hand. 

Titan  Research  and  Technology's  (TRT)1  CRALE  1-  and  2-D  finite  difference 
codes  assume  the  stress  tensor  is  separable  into  the  mean  stress,  P,  and  a  deviatoric 
stress  tensor,  ay'.  The  algorithms  used  to  obtain  the  mean  stress  are  collected  in 
subroutines  commonly  referred  to  as  the  Equation  of  State  (EOS).  The  separation  of  P 
and  ay'  allows  the  code  to  bypass  calculations  of  stresses  and  strains  for  those 
,  materials  that  can  only  support  a  pressure,  e.g.,  gases  and  fluids  such  as  air,  water  or 
high  explosive  (HE),  or  a  solid  which  has  melted.  Such  materials  only  require  the 
current  density  and  energy  to  obtain  the  pressure.  For  more  complex  material 
behavior,  the  code  assumes  that  the  stress  deviators  can  be  calculated  from  the 
previous  stresses  and  current  increments  of  the  strain  tensor,  Aey  and  that,  at  most, 
only  four2  deviatoric  stresses  are  required,  i.e.,  a'xx,  o'yy,  a'^,  a'Xy.  Finally,  provisions 
are  included  to  account  for  hysteresis,  fracture,  and  other  nonelastic  responses  such 
as  solid-solid  phase  changes,  plastic  flow  and  dilation. 

The  most  general  constitutive  relations  in  CRALE  model  the  behavior  of  an 
isotropic,  elastic-plastic,  hysteretic  material  which  can  crack  under  tensile  loads.  For 
simpler,  purely  hydrodynamic  models,  e.g.,  HE  and  WATER,  no  deviatoric  stresses  are 

^  A  division  of  Titan  Corp,  formerly  California  Research  and  Technology,  CRT 

2  The  missing  stress  components  are  either  zero  due  to  the  2D  symetry  or  equal  to  o'yy  due  to  equilibrium 
requirements. 
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present  and  the  material  pressure  state  is  completely  specified  by  casting  the  pressure 
as  a  function  of  the  current  density  (p.RHO3)  and  specific  energy  (E ,ENG).  [Rather 
than  density,  CRALE  passes  the  excess  compression,  p  (=  [p/p0]-1 ,  EMUL)  to  minimize 
roundoff  error.] 

Since  a  wide  range  of  materials  and  loading  conditions  may  exist  in  any  specific 
calculation,  a  number  of  different  EOS  routines  are  available.  CRALE1  and  2  have 
been  written  so  that  general  routines  (SHEL,  MIX,  VIC,  CIST,  MCIST,  HE,  WATER,  and 
CAP)  are  available  to  model  a  wide  variety  of  materials.  In  addition,  the  EOS 
calculation  is  linked  to  the  main  codes  by  passing  variables  through  named  common 
blocks,  /EOS/,  /EOS2/  and  /EOS3/,  so  that  the  user  may  include  his  own  routine  using 
a  shell  routine,  UEOS,  with  a  minimum  of  effort.  (Additional  EOS  models  have  been 
incorporated  on  various  program  libraries,  PL's,  from  time  to  time,  e.g.  AFWL,  S3EOS, 
WAGNER  and  GRAY.  These  routines  have  been  removed  for  the  current  codes  as 
obsolete,  but  could  easily  be  reactivated  if  required.) 

1.2  REPORT  STRUCTURE. 

This  report  is  intended  to  provide  a  description  of  the  routines  used  by  the 
CRALE  codes  to  obtain  the  stress  state  in  a  cell.  A  glossary  of  terms  used  in  the  codes 
follows  in  Section  1 .3.  A  general  overview  of  the  material  modeling  algorithms  used  in 
CRALE  is  presented  in  Section  2.  A  summary  of  the  general  input  required  by  the 
EOS  models  and  the  various  EOS  models  currently  available  in  CRALE  are  described 
,  in  Section  3. 

\ 

As  part  of  the  procedures  used  during  development  of  the  material  models,  a 
program,  EQS,  was  written  that  allows  the  user  to  drive  any  EOS  module  through  a 
variety  of  stress-strain  paths  and  compare  the  results  with  experimental  data.  A  User's 
guide  for  EQS,  with  detailed  descriptions  of  the  I/O  options,  is  presented  in  Section  4. 

Various  data  bases  have  been  constructed  and  maintained  for  DNA  by  TRT  for 
use  in  constructing  and  verifying  the  EOS  models  and  to  provide  archival  storage  for 
data.  These  data  bases  are  discussed  in  Section  5.  Instructions  to  allow  the  DNA  user 
community  access  to  these  files  is  also  included. 


I 

-1 _ 

®  The  code  names  for  the  various  quantities  are  shown  in  red  and  Italic.. 
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1.3  SYMBOLS  AND  DEFINITIONS. 

The  fundamental  units  in  the  CRALE  codes  and  the  EOS  routines  are 
gram-centimeter-microsecond  (g-cm-|xs).  Hence,  the  units  of  pressure  (stress)  and 
energy  are  Megabars  (Mb)  and  Terraergs  (Te),  respectively.  A  glossary  of  terms  used 
in  this  report  is  presented  in  Table  1-1 .  Model  specific  variables  are  defined  in  Section 

3. 


Table  1  -1 .  Glossary  of  general  EOS  terms  and  symbols. 


Variable  Units  Code  Symbol  Definition 


B,  K 

Bo 

Bm 

Cl 

Cp 

Cs 

E 

Em-Bs 

G 

I 

r 


Mb 

Mb 

Mb 

cm/|is 

cm/ps 

cm/|is 

Te/g 

Te/g 

Mb 


BO 

BM 


ENG 

EM.ES 

G 


M 

>2 

*3 


J 

Mb 

J' 

Mb 

Jl 

Mb 

J2 

Mb2 

J'2 

Mb2 

J2P 

^3 

Mb3 

P 

Mb 

P 

Y 

Mb 

Yo 

Mb 

YO 

Yvm 

Mb 

YVM 

p 

Te/g 

CTE 

r 


^  EMU 

e 


£ii 

£,ii 


v 


p 

g/cm3 

Po 

g/cw? 

RHOZ 

Pr 

g/cm^ 

RHOREF 

a 

Mb 

°ii 

Mb 

°'ii 

Mb 

Bulk  moduli,  dP/dp 
Initial  Bulk  moduli 
Maximum  Bulk  moduli 

Longitudinal  sound  speed,  Cl2  =  (K  +  4/3G)/po 

Bulk  sound  speed,  C02  =  K/p0 

Shear  wave  speed,  Cs2  =  G/p0 

Specific  internal  energy  of  the  material 

Specific  internal  energy  of  vaporization  of  the  material 

Shear  modulus 

The  strain  tensor;  i.e.,  ejj,  i,j=1,3. 

The  deviatoric  strain  tensor;  i.e.,  e'jj,  i,  j=1l3. 

1st  invariant  of  the  strain  tensor  =  Eeji,  i=1,3 
2nd  invariant  of  the  strain  tensor;  =  1/2Lejjejj,  i,j=1,3 
3rd  invariant  of  the  strain  tensor;  =  |ejj|,  i„j-1.3 
The  stress  tensor;  i.e.,  Ojj,  i„j=1,3. 

The  deviatoric  stress  tensor;  i.e.,  a'jj,  itlj=1,3. 

1st  invariant  of  the  stress  tensor;  =  3P  =  lojj,  i®1,3 
2nd  invariant  of  the  stress  tensor;  =  1/2£ojjajj,  i,  js1,3 
2nd  invariant  of  the  deviatoric  stress  tensor;  =  I^Eo’jjc’jj,  i,,j=1»3 
3rd  invariant  of  the  stress  tensor;  =  |<*jj|; 

Pressure  or  mean  stress,  =  Ji/3 

Plastic  yield  surface,  =  the  maximum  elastic  (J'2 )1/2 

Cohesion,  the  yield  surface  at  P=0 

Maximum  Y;  named  for  Von  Mises 

Coefficient  of  thermal  expansion 

Y-1  in  perfect  gas  EOS 

Compression;  =  p/po 

Excess  compression;  =  t|~1  =  p/po*1 

Strain 

The  ij  component  of  the  strain  tensor 

The  ij  component  of  the  deviatoric  strain  tensor;  =  z\\  -  J^5jj/3 

Poissons  Ratio 

Density  of  the  material 

Initial  density  of  the  material 

Reference  density  of  the  material 

Stress 

The  ij  component  of  the  stress  tensor 

The  ij  component  of  the  deviatoric  stress  tensor;  =  ojj  -  PSjj 
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During  the  course  of  writing  this  report,  it  became  obvious  that  a  number  of 
terms  in  common  usage  appear  to  be  somewhat  loosely  defined.  For  the  purposes  of 
this  report  the  following  conventions  are  assumed. 

•  Stresses  and  strains  are  considered  positive  in  compression.  Note,  this 
convention  is  consistent  with  P  and  p,  but  generally  the  opposite  of 
standard  engineering  usage. 

•  The  term  "Equation-of-State,  EOS"  is  used  in  the  narrow  sense  of  relating 
the  mean  stress  (pressure)  to  the  density  and  energy.  "Constitutive 
model"  refers  to  the  complete  stress-strain  relation.  Unfortunately,  the 
terms  are  frequently  used  interchangeably.  EOS  was  commonly  used  by 
physicists  to  describe  high  temperature,  high  pressure  states;  engineers 
interested  in  the  low  stress  behavior  of  materials  tended  to  talk  of 
constitutive  behavior.  Since  calculations  in  the  codes  frequently  require 
models  that  represent  a  continuum  of  behavior  spanning  pressures 
ranging  from  tensions  to  multi-Megabars,  the  distinctions  between 
Equations  of  State  and  constitutive  models  have  become  blurred. 

•  Similarly,  the  terms,  elastic,  non-elastic  and  plastic  need  to  be  defined. 

We  use  elastic  to  describe  incremental  changes  in  deviators  that  obey 
Hook's  law,  i.e.,  6a  =  M8e,  where  the  modulus,  M,  is  only  linear  in  the 
increment.  Any  other  changes  to  the  stress  tensor  are  non-elastic,  e.g., 
fracture,  hysteresis  or  phase  changes.  Plasticity  describes  only  those 
non-elastic  changes  that  result  when  deviatoric  stresses  try  to  exceed  the 

.  strength  of  the  material. 

•  Strains  are  incrementally  computed  from  displacement  and  are  the  natural 
strains  (Av/v),  not  small  strain  or  engineering  strain.  Small  strain  ignores 
the  cross-products  in  computing  the  volumetric  strain;  engineering  strain 
relates  changes  in  displacement  to  the  initial  rather  than  current  positions. 
Differences  between  these  definitions  become  significant  at  strains  of 
about  5%.  Transformations  between  the  strain  definitions  are  possible.  It 
should  be  noted  that,  due  to  non-elastic  effects,  the  differences  between 
the  various  strains  do  not  translate  into  similar  changes  in  stress.  In  fact, 
small-strain  codes  have  frequently  generated  accurate  solutions  even 
when  material  underwent  extreme  distortion. 

•  The  terms  'pressure'  and  'mean  stress'  are  used  interchangeably. 


i 
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SECTION  2 


GENERAL  EQUATION  OF  STATE  ALGORITHMS 

CRALE,  typical  of  finite  difference,  continuum  mechanics  codes,  requires  values 
of  the  complete  constitutive  behavior  of  the  material  in  each  computational  cell  for  each 
cycle  in  order  to  derive  a  solution.  The  EOS  models  which  relate  the  thermodynamic 
variables,  i.e.,  stress  and  energy,  to  the  mechanical  variables,  position  and  density  and 
to  the  past  history  of  the  material  provide  one  part  of  the  behavior,  additional  algorithms 
are  needed  to  complete  the  description. 

In  general,  the  calculation  of  the  constitutive  behavior  for  each  cell  in  a  CRALE 

solution  proceeds  in  5  steps,  namely; 

1 .  update  the  cell  coordinates  based  on  previous  positions  and  velocities, 

2.  calculate  the  density,  p,  and  specific  energy,  E, 

3.  call  the  appropriate  EOS  module  to  obtain  the  new  mean  stress  and 
elastic  moduli,  K  and  G 

4.  if  necessary,  compute  incremental  changes  in  the  strain  tensor,  Ae'y, 

and  calculate  an  intermediate  stress  tensor  assuming  the  changes  are  elastic, 

5.  test  and  adjust  the  stress  tensor,  as  needed,  based  on  plastic  yielding, 
tensile  failure. 

A  more  complete  description  of  the  general  methodology  of  steps  3  through  5 
'follows.  Model  specific  differences  from  these  algorithms  are  presented  in  the 
'discussions  of  the  various  EOS  models  in  Section  3. 

2.1  PRESSURE. 

In  general,  CRALE  assumes  that  the  mean  and  deviatoric  components  of  stress 
are  separable.  Furthermore,  the  mean  stress  is  composed  of  two  parts;  a  solid  and  a 
vapor  component.  Both  components  are  functions  of  p  and  E,  but  the  solid  portion  is 
primarily  a  function  of  density  while  the  vapor  pressure  is  governed  by  the  energy  per 
unit  volume  (pE).  The  CRALE  code  provides  the  various  EOS  models  with  the  current 
density,  internal  energy  and  a  history  parameter  and  expects  the  pressure,  sound 
speed,  bulk  modulus  and  shear  modulus  in  return.  Any  coupling  of  the  pressure  to  the 
deviatoric  stresses  is  included  in  the  history  parameter.  Currently,  the  maximum 
Excess  compression,  (j.max,  which  the  material  has  experienced  is  used  as  the  history 
parameter  primarily  to  determine  the  effect  of  hysteretic  compaction  on  the  pressure 
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and  moduli.  To  minimize  potential  round-off  errors,  p,  the  excess  compression  is 
carried  by  the  codes,  rather  than  p. 

2.2  ELASTIC  DEVIATORIC  BEHAVIOR. 

For  materials  that  support  deviatoric  stresses,  the  EOS  module  returns  a  non¬ 
zero  shear  modulus  to  the  main  code.  The  code  then  calculates  the  change  in  the 
strain  tensor  and  computes  the  incremental  elastic  change  in  the  deviatoric  stresses  by 
the  simple  Hook's  law  Equation; 

AOy  =  -2GAeJj  (2.1) 

The  negative  sign  is  required  because,  in  CRALE,  strains  are  defined  as  positive 
in  tension  while  stresses  are  positive  in  compression.  The  shear  modulus  is  calculated 
in  the  EOS  each  cycle  and  may  vary  with  compression  or  pressure  so  that  the  stresses, 
in  general,  are  only  incrementally,  linearly  elastic. 

2.3  PLASTIC  YIELDING. 


The  CRALE  codes  assume  that  the  deviatoric  stresses  a  solid  material  is 
capable  of  supporting  are  limited  by  a  surface  in  stress  space  which  is  a  function  of  the 
mean  stress  and  the  sign  of  the  third  invariant  of  the  deviatoric  stress  tensor,  J3 

although  more  sophisticated  surfaces  such  as  that  incorporated  in  the  CAP4  model  are 
also  possible.  The  surface  currently  used  in  the  codes  is  defined  by  relating  the 
,  maximum  allowable  second  invariant  of  the  deviatoric  stress  tensor,  J2,  and  the 

^  pressure,  thus; 

V^<F(P)  =  F(J,/3)  (2.2) 

Yvm-(Yvm-Y0)e'<P,Py) 


F(P)  = 


or 


max  [Yvm,Y0  +  YjP] 


(2.3) 


where  Yvm,  Y0  and  Pi  are  input  parameters  which  may  depend  on  J3  and  whether  the 
material  is  intact  or  fractured. 


I 

J _ 

4  G.Y.  Baladi,  "An  Effective  Stress  Model  for  Ground  Motion  Calculations",  WES  TR  SL-79-7,  Sept  1979 
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In  CRALE2,  the  calculations  of  the  deviatoric  stress  tensor  are  performed  in 
subroutine  STRAIN,  so  the  EOS  subroutines  only  set  the  elastic  moduli  and  elastic 
yield  surface  parameters  for  the  current  material.  The  yield  parameters  are  stored  in 
CK(28,NE)  through  CK(37,NE)  for  all  the  EOS  models  and  are  illustrated  in  the  sketch 
below.  As  shown  in  the  figure,  the  yield  surface  is  generally  of  the  Drucker-Prager  form 
with  a  von  Mises  limit,  although  strain  hardening  or  softening  of  the  initial  yield  surface 
may  be  included  by  setting  the  hardening  and  softening  parameters  ( BETAH  and 
BETAS)  and  the  value  of  plastic  strain  ( EPYLD )  at  which  the  material  stops  hardening 
and  begins  to  soften.  The  yield  surface  can  never  harden  above  YLDVM  or  soften 
below  the  surface  for  3-cracked  material. 


If  exceeds  the  yield  surface  the  code  imposes  a  flow  rule  that  adjusts  the 
deviators  so  that  equals  Y.  Two  basic  options  exist;  associated  and  non- 
associated  flow  .  The  non-associated  flow  simply  reduces  the  stress  deviators 
proportionally  without  changing  the  mean  stress,  i.e.,  along  a  path  perpindicular  to  the 
pressure  axis  in  the  sketch  and  is  the  default  condition.  Associated  flow  (activated  by 
, setting  YFLOWR="\)  moves  4^2.  along  a  line  perdindicular  to  the  yield  surface  itself. 
This  adjustment  produces  a  change  in  the  mean  stress  along  with  the  change  to  the 
deviators  and  can  generate  nonphysical  states. 

2.4  FRACTURE. 

The  CRALE  code  allows  material  to  crack  under  several  conditions.  After 
completing  the  calculation  of  the  pressure  and  stress  deviators  assuming  the  material 
is  intact,  the  stress  tensor  is  rotated  to  obtain  the  three  principal  stresses,  Ojj  j=-j  3. 

The  code  allows  cracks  to  open  in  several  possible  ways,  i.e.,  as  a  single  crack,  as  two 
orthogonal  in-plane  cracks,  as  a  combination  of  an  in-plane  and  a  hoop  crack  or  in  all 
directions. 

,  To  determine  if  a  previously  intact  cell  is  cracked,  each  Cj  is  compared  to  the 
tensile  limit,  cr  If  any  one  principal  stress  is  less  than  at  (at<0),  a  crack  is  opened. 
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The  offending  stress  is  set  to  zero  and  the  other  two  stresses  are  adjusted  assuming 
the  change  is  due  to  a  change  in  the  strain  perpendicular  to  the  Oj  direction.  The 
principal  angle,  ©j(  is  saved  and  the  remaining  two  stresses  rechecked  to  insure  that 
they  are  still  greater  than  ot.  The  amount  of  strain  required  to  return  the  stress  to  zero 
is  also  stored  as  a  crack  volume  ( VCR1 ).  If  the  stress  adjustment  has  driven  a  second 
principal  stress  beyond  ct,  the  code  also  sets  it  to  zero  and  resets  the  third  stress, 
again  assuming  a  uniaxial  strain  change,  and  saving  the  adjustment  as  a  second  crack 
volume  ( VCR2 ).  On  subsequent  cycles,  the  material  volume  is  reduced  by  the  sum  of 
VCR1  and  VCR2  when  computing  the  pressure. 

After  a  cell  has  cracked,  the  code  rotates  the  stresses  using  ©j  as  one  reference 
angle;  no  shear  stress  is  allowed  parallel  to  this  direction.  Once  one  or  two  cracks 
have  opened,  the  stress  perpendicular  to  the  crack  is  continually  returned  to  zero  and 
the  crack  volume  adjusted.  Upon  recompression,  the  stress  is  not  permitted  to  go 
positive  until  the  crack  volume  has  been  reduced  to  zero. 

If  all  three  principal  stresses  exceed  the  tensile  limit,  the  material  is  assumed  to 
have  shattered.  In  this  case,  all  the  stresses  and  the  crack  volumes  are  set  to  zero, 
hence  the  material  can  only  support  positive  stresses  upon  recompression  if  the 
pressure  is  positive. 

Because  of  its  unique  treatment  of  stress,  the  fracture  logic  is  bypassed  for 
material  using  the  CAP  model  EOS. 
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SECTION  3 


EOS  MODELS 

3.1  SUBROUTINE  EQST:  OVERVIEW  AND  GENERIC  INPUT. 

In  order  to  permit  some  flexibility  in  using  multiple  material  models  in  a  single 
calculation,  all  of  the  EOS  routines  are  accessed  by  CRALE1 ,  CRALE2  and  their 
ancillary  programs  via  calls  to  one  main  subroutine,  EQST.  Hence  adding  new  models 
or  changing  the  EOS  I/O  only  requires  updating  EQST.  The  material  properties  for  the 
models  needed  are  initialized  by  a  call  to  ESINIT,  an  entry  point  in  EQST.  After  all  the 
model  data  have  been  read  in,  the  parameters  are  echoed  to  the  output  file  by 
ESPRNT,  a  second  entry  point  in  EQST. 

The  parameters  of  the  various  EOSs  are  obtained  by  the  code  in  several  ways. 
In  general,  the  code  initializes  the  EOS  parameters  via  a  CALL  to  ESINIT,  an  entry 
point  in  subroutine  EQST.  ESINIT  reads  an  initial  line  with  the  name  of  the  material, 
three  flags  and  a  48-character  comment.  The  first  flag,  NETYPE,  designates  which 
EOS  model  to  use;  the  second,  NEO,  if  the  input  data  is  to  be  read  from  the  data  base 
on  Tape  14;  and  the  third,  ICNES,  the  units  of  the  EOS  input  parameters  for  this 
material.  If  the  material  is  not  available  from  the  data  base,  i.e.,  NEO= 0,  the  code 
proceeds  to  read  a  second  line  with  the  initial  ( RHOZ)  and  reference  ()  densities,  the 
minimum  vaporization  energy  density  (EM),  the  initial  energy  density  (EZ),  and  the 
,  tensile  limit- {ST).  Knowing  which  EOS  to  use,  the  code  now  reads  the  necesary 
additional  parameters  mto  the  CK  block.  The  number  of  lines  of  input  parameters 
varies  with  each  EOS,  the  more  complicated  models  requiring  more  data.  Currently, 
subroutines  SHEL,  MIX,  CIST,  and  CAP  read  40  input  parameters  (i.e.  5  data  line  with 
8  numbers  per  line);  HE  only  reads  two  additional  lines,  16  variables.  The  parameters 
required  by  each  EOS  are  discussed  in  Section  4.  As  shown  schematically  in  Table 
3-1  and  defined  in  Table  3-2,  data  sets  for  up  to  15  different  materials  can  be  read,  in 
arbitrary  order,  for  any  single  code  run.  The  input  file  for  phyllite  (as  echoed  by 
ESPRNT)  is  included  in  Table  3-1  as  an  example  for  a  material  using  the  SHEL  EOS 
model.  Similar  input  files  can  be  recalled  from  the  EOSS.dat  file;  a  list  of  which  is 
provided  in  Table  5-8.  This  data  base  can  be  accessed  by  the  codes  for  ease  of 
inputing  specific  materials  and  to  insure  uniformity  throughout  the  various  programs. 
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All  the  EOS  used  by  the  CRALE  codes  and  their  auxiliary  routines  use  the 
named  COMMON  blocks,  /EOS1/,  /EOS2/  and  /EOSY/  to  store  and  transfer  the  values 
of  the  parameters  and  the  input/output  quantities  required  for  each  call  to  the  EOS. 
Table  3-3  lists  each  of  these  variables  with  its  size,  units,  and  definition.  An  I  or  0  in 
the  column  labeled  I/O  signifies  this  variable  is  either  passed  into  (I)  or  calculated  and 
returned  from  (0)  the  EOS  every  call.  The  other  quantities  in  the  COMMON  blocks 
remain  constant  for  each  material  or  for  the  entire  problem  (e.g.,  ATMOS). 

Table  3-1 .  Summary  of  input  for  the  EOS  models. 


LINE 

VARIABLES  / 

FORMATS 

1 

IEOS 

A6 

NETYPE 

12 

NEO 

12 

ICNES 

12 

TMATRX 

A6 

TSOLID  TFLUID  WORDS 

A6  A6  12A4 

la 

CONV 

8F10.0 

(if  needed) 

1b 

UNITS 

8A8 

(if  needed 

2 

RHOZ 

F10.0 

RHOREF 

F10.0 

ES 

F10.0 

EZ 

F10.0 

ST 

F10.0 

3 

CK1-N 

8F10.0 

N  =40  for  SHEL,  CIST,  MIX,  CAP 
=16  for  HE 

H 

END 

LIKE  1 
LIKE  2 
LIKE  3a 
LIKE  3b 
LIKE  3c 
LIKE  3d 
LIKE  3e 


PHYLL  1  an  example  of  the  input  for  the  SHEL  EOS 
2 . 7 800E+00  2 . 7  800E+00  7.0000E-02  O.OOOOE+OO  6.6700E-05 
4. 0000E-01  8 . 0000E-01  2.0000E-01  O.OOOOE+OO  O.OOOOE+OO  O.OOOOE+OO 
0 . OOOOE+OO  O.OOOOE+OO  O.OOOOE+OO  O.OOOOE+OO  O.OOOOE+OO  O.OOOOE+OO 
3 . OOOOE+OO  O.OOOOE+OO  O.OOOOE+OO  4.2500E+00  2.0000E-02  1.3000E-01 
3 . OOOOE-Ol  O.OOOOE+OO  O.OOOOE+OO  O.OOOOE+OO  O.OOOOE+OO  O.OOOOE+OO 
8 . 0000E-03  5 . 0000E-04  2.5000E+02  O.OOOOE+OO  O.OOOOE+OO  O.OOOOE+OO 


O.OOOOE+OO 
7.5000E-02 
3 . 5000E+00 
O.OOOOE+OO 
O.OOOOE+OO 


O.OOOOE+OO 
O.OOOOE+OO 
2. OOOOE+OO 
O.OOOOE+OO 
O.OOOOE+OO 
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Table  3-2.  EOS  input  variables. 


VARIABLE 


LME=  1  FORMAT  =  A6,  I2,2x,I2,3A6,12A4 

IEOS  Name  of  the  Equation-of-State  (EOS) 

NETYPE  Number  of  the  EOS  Subroutine 
1  -  SHEL  (CRT) 

2 -VIC  (CRT) 

3  -  CIST  (  Modified  AFWL  EOS  ) 

4  -  HE  (  LLL  Hi  Explosive  ) 

5  -  MCST  (CRT) 

6  -  SCUB  (  Obs  SCUBED  Table  Lookup  ) 

7  -  H20  Walker-Srernberg  Water 

8  -  WAGN  (Obs  CRT  model) 

9  -  CAP  (WES  form  of  the  CAP  EOS) 

10- MIX  (CRT) 

NEO  if  >  0  Read  EOS  input  for  this  material  from  Tape  14 
ICNES  if  =  0  EOS  data  in  standard  code  units 

=  1-8  use  ICNES  set  of  units  from  table  for  this  EOS  input 
=  9-10  read  in  (lines  1  a  &  1  b)  and  use  special  set  of  units  for  this  input 
TMATRX  the  matrix  EOS  used  if  INEOS=1 0 
TSOLID  the  solid  EOS  used  if  INEOS=10 
TFLUID  the  fluid  EOS  used  if  INEOS=10 
WORDS  Comment 


LIME  =  la  FORMAT  =  8F10.0  Optional:  read  only  if  ICNES  >  8 

CONV  New  set  of  conversion  factors  for  EOS 
LINE  =  1b  FORMAT  =  8A8  Optional:  read  only  if  ICNES  >  8 

UNITS  Names  of  new  units 


LINE  =  2  FORMAT  =  8F10.0 

RHOZ  Initial  density  (p0)  of  material 

RHOREF  Reference  density  (pref)of  material 

ES  Max  energy  of  solid  material 

EZ  Initial  energy  of  zones 

ST  Tensile  limit,  <0 


LINES  =  3  FORMAT  =  8F1 0.0 

CK  1-N  Parameters  of  the  EOS.  N  varies  according  to  the  type  of  EOS. 

See  Sec  4  for  details  of  EOS  input. 

UME=4  FORMAT  =  A6 

i  END  EOS  data  terminated  by  line  with  IEOS  =  3hEND 
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Table  3-3.  General  parameters  used  by  all  Equations-of-state. 


VARIABLE 

UNITS 

DESCRIPTION 

ATMOS 

Mbar 

Atmospheric  pressure,  generally  set  to  1  bar 

BETAH 

Mbar 

Strain  hardening  parameter 

BETAS 

Mbar 

Strain  softening  parameter 

CK 

Parameters  used  by  the  EOSs 

COHC 

Mbar 

material  cohesion  in  compression 

COHE 

Mbar 

material  cohesion  in  extension 

COHES 

COHF 

Mbar 

cohesion  of  3-cracked  material 

CSS 

cm/s2 

Sound  speed  of  cell 

DETA 

— 

The  change  in  EMUL  on  this  cycle 

DMUU 

— 

Maximum  p  of  cell  (pmax) 

DPDE 

Mb-g/Te 

Derivative  of  P  with  energy 

DPDMU 

Mbar 

Derivative  of  P  with  compression 

DYDPC 

DY/DP 

of  material  in  compression 

DYDPE 

DY/DP 

of  material  in  extension 

DYDPF 

DY/DP 

of  3-cracked  material 

EFRAC 

Te/g 

E  -  Em)/Em 

EM 

Te/g 

Set  to  ES(NE) 

EMUL 

— 

Excess  compression  of  cell 

ENG 

Te/g 

Specific  internal  energy  of  the  cell 

EPSKK 

EQUATIOND 

UM 

EPYLD 

Plastic  strain  where  Y  switches  from 

ES 

Te/g 

strain  hardening  to  softening 

Vaporization  energies  of  each  EOS 

EZ 

Te/g 

Initial  energy  densities  of  each  EOS 

G 

Mbar 

Shear  modulus 

'  GEOP 

Mbar 

Geostatic  pressure  in  cell 

ICNES 

Allows  input  in  units  other  than  standard 

'  ICR 

— 

Crack  state 

IEOS 

— 

Debug  sentinel  (with  ‘DEFINE  DEBUG) 

IEOSPR 

— 

NE 

— 

Number  of  the  EOS  for  this  material 

P 

Mbar 

Pressure 

PAIR 

Mbar 

The  air  pore  pressure 

PM 

Mbar 

Pressure  on  last  cycle 

PYLD 

Mbar 

Obsolete 

RHOREF 

g/cm3 

Reference  densities  for  each  EOS 

RHOZ 

g/cm3 

Initial  densities  for  each  EOS 

ST 

Mbar 

Tensile  strengths  for  each  EOS 

STRSS 

— 

Reference  specific  vol  of  each  EOS 

VZ 

cm3/g 

YFLOWR 

Flow  rule  sent  (0=non-assoc,1=assoc) 

YLDMXC 

Mbar 

Maximum  Y  in  compression 

1  YLDMXE 

Mbar 

Maximum  Y  in  extension 

1  YLDMXF 

Mbar 

Maximum  Y  for  3-cracked  material 
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3.2.  SPECIFIC  EOS  MODELS. 


The  various  subroutines  described  in  this  section  generally  calculate  the  mean 
stress  or  pressure  (P),  bulk  and  shear  moduli  (K,G)  and  longitudinal  sound  speed 
(CSS)  for  any  pair  of  input  values  of  density  (p),  specific  energy  (E)  and,  for  hysteretic 
materials,  the  maximum  compression  (pmax).  Elastic  moduli  are  only  required  for 
materials  that  support  deviatoric  stresses  and  can  be  set  to  zero  for  gases  or  fluids.  (If 
G  is  set  to  zero  the  code  will  suppress  the  calculation  of  deviatoric  stresses  and  strains 


The  parameters  used  by  the  various  EOS  subroutines  are  stored,  by  material,  in 
the  two-dimensional  array  CK(i,NE),  summarized  in  Table  3-4,  and  defined  in  the  tables 
accompanying  the  description  of  each  model  that  follows.  For  each  of  the  CK  variables 
available  to  the  specific  EOSs,  the  mnemonic  eqivalent  is  given,  followed  by  a  Y  (yes) 
in  the  INPUT  column,  if  this  CK  is  an  input  number.  If  the  CK  is  modified  by  the 
subroutine  after  input  and  replaced  by  a  new  value  stored  in  the  same  location,  a  Y 
appears  in  the  MOD  (modified)  column;  a  C  in  this  column  indicates  the  value  is 
derived  from  other  input.  The  internal  code  units  of  the  parameter  are  given  next, 
followed  by  a  short  description.  Parameters  calculated  for  information  only  (FIO)  are  so 
labeled.  These  may  include  data  which  are  no  longer  used  by  the  code  (such  as 
PYLD)  and  derived  numbers  which  may  be  useful  (e.g.,  initial  wave  speed  and 
porosity).  Only  the  first  50  CK  variables  are  echoed  to  the  output  file. 

The  SHEL  EOS  model  was  originally  developed  in  the  late  1960's  when  the 
'  scope  of  many  finite  difference  code  solutions  was  expanded  to  include  pressures  and 
\  stresses  that  extended  from  tenths  of  a  bar  to  many  Megabars  with  particular  emphasis 
on  ground  media,  i.e.,  soils  and  rocks.  A  continuos  EOS  model  capable  of  including 
elastic  and  non-elastic  behavior  as  well  as  the  high  temperature  pressure-volume- 
energy  relations  was  needed.  Over  the  years,  the  SHEL  model  has  evolved  to  include 
several  hysteretic  options,  solid-solid  phases  and  other  changes  needed  to  describe 
specific  material  behavior.  Several  of  the  other  models,  i.e.,  VIC,  MCST,  MIX  use  or 
are  derivatives  of  the  SHEL  EOS.  Hence,  the  SHEL  EOS  routine  is  discussed  in  detail 
in  Section  3.2.1;  the  other  model  descriptions  reference  it  as  appropriate. 


J 


l 
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Table  3-4.  Summary  of  CK  block  parameters  of  current  EOS  models. 


EOS 

# 

NAME 

CKs 

i 

SHEL 

2 

VIC 

3 

CIST 

1 

BO 

CUR 

Cl 

2 

BM 

Cl 

C2 

3 

XMUO 

C3 

4 

PTOE 

PI 

PI 

5 

A1 

P2 

6 

BTOE 

P3 

P3 

7 

AMU1 

XMU 

EPOT 

8 

A2 

ALFA 

9 

DMUMN 

EMU3 

EMU3 

10 

DMUMX 

ALPH1 

EMU2 

11 

EMU1 

EMU1 

12 

EMUMX 

EMUMAX 

13 

B1ELAS 

14 

BUR 

B1 

15 

STOV 

DPDMU1 

B2 

16 

SWTH 

DPDMUK 

B3 

17 

CTE 

CTE 

CTE 

18 

ENGPH 

PMU2 

EV 

19 

PHI 

FPRL 

FPRL 

20 

ZMU2 

FPRU 

FPRU 

21 

PHRE1 

ALLPHA 

22 

PHRE2 

PMU1 

23 

BM2 

PMU2 

24 

PHC 

BETA 

25 

POIS 

PRL 

26 

FPOIS 

PRU 

27 

GMAX 

GMAX 

GMAX 

28 

BETAH 

BETAH 

BETAH 

29 

BETAS 

BETAS 

BETAS 

30 

EPYLD 

EPYLD 

EPYLD 

31 

YFLOWR 

YFLOWR 

YFLOWR 

32 

YLDMXF 

YLDMXF 

YLDMXF 

33 

YLDMXC 

YLDMXC 

YLDMXC 

34 

COHC 

COHC 

COHC 

35 

DYDPC 

DYDPC 

DYDPC 

36 

YLDMXE 

YLDMXE 

YLDMXE 

37 

COHE 

COHE 

COHE 

38 

DYDPE 

DYDPE 

DYDPE 

39 

COHF 

COHF 

COHF 

40 

DYDPF 

DYDPF 

DYDPF 

41 

EMUZ 

42 

ETAZ 

43 

EMU1 

44 

45 

46 

47 

48 

49 

VREF 

VREF 

VREF 

50 

NEOS 

NEOS 

NEOS 

4 

5 

8 

9 

HE 

MCST 

UEOS 

CAP 

HEGAM 

A 

user 

AKEI 

A 

B 

specified 

AK1 

R1 

EO 

model 

AK2 

B 

ESS 

AK3 

R2 

ESP 

AK4 

DETV 

ALPHA 

AK5 

AGAMM 

BSQ 

AK6 

CJP 

AA 

AK7 

AGAM 

BB 

AKIM 

GAMPG 

GM 

AKIM 

POIS 

AK2M 

GOK 

AA 

DELZ 

AB 

AB1 

AC 

A1 

AP 

AGE! 

A2 

BP 

AG1 

CUTA 

EPM 

AG2 

B1 

EPP 

AG3 

B2 

DEP 

AG4 

CUTB 

CRO 

CUTC 

CR1 

CR2 

AW 

AD 

ADI 

BETAH 

AD2 

BETAS 

AD3 

EPYLD 

AD4 

CJMU 

YFLOWR 

CJGAM 

YLDMXF 

YLDMXC 

SJMX 

COHC 

CKMX 

DYDPC 

CKAA 

YLDMXE 

LTYPE 

COHE 

DYDPE 

COHF 

DYDPF 

ASTART 

XSTART 

FCUT 

TCUT 

VREF  VREF  VREF  VREF 

NEOS  NEOS  NEOS  NEOS 


I 


I 


10 

MIX 


c 

B 

GAMZ 

AGAM 

GAMMZ 

GAME 

GP 

TMZ 

VJ 

VB 

AYP 

DELSP 

THETA 

ALFAP 

RP3 

PCC 

EZZ 

DSA 

SMA 

ZJ 

XJ 

RP3SQ 


BETAH 

BETAS 

EPYLD 

YFLOWR 

YLDMXF 

YLDMXC 

COHC 

DYDPC 

YLDMXE 

COHE 

DYDPE 

COHF 

DYDPF 


VREF 

NEOS 


14 


3.2.1  TheSHELEOS. 

3.2. 1.1  General. 


Given:  a  state  of  the  material  defined  by  its  current  p,  ey;  E  and  some  measure  of  the 
past  strain  path 

Find:  the  pressure  or  mean  stress,  P  and  the  deviatoric  stress  tensor,  a  y 
Definitions  for  the  parameters  used  in  SHEL  are  presented  in  Table  3-5. 

The  SHEL  EOS  proceeds  as  follows.  Assume  the  pressure  is  composed  of  two 
terms,  a  solid  part,  Ps,  which  depends  on  compression,  energy  and  a  past  history 
parameter  and  a  vapor  part,  Py,  which  is  a  function  of  p  and  E.  Then, 


P(p,E)  =  Ps(p,E)+Pv(p,E)  V-V 

Ps  =  B^MBm  -B0)n*  (l-  e-^')  <32> 

Pv  =  rpE*;  (E>Em)  (33) 

where 

p'=p  +  pE  <3-4) 

E*  =  (E  -  Em  )(1  -  e‘tE'Em]/Em )  <3-5) 

sp  =  y  - 1  =  0. 6  +  0. 023[ln(p  / E)f  +  0. 1 4  •  ln(p  /  E) + 0. 05  •  In  p  (3.6) 


Five  parameters  are  required  to  compute  the  pressure  (Equation  3.1)  in  a 
material  with  no  air-porosity  (voids),  i.e.  p0=pr,  'n  addition  to  the  reference/initial  density. 
'  They  are: 

Em  =  minimum  energy  required  to  activate  Pv,  (Te/g) 

B0  =  the  zero  pressure  bulk  modulus,  dP/dp  (Mb) 

Bm  =  the  maximum  bulk  modulus,  (Mb) 
p*  =  the  excess  compression  at  which  the  bulk  modulus 
has  increased  exponentially  to  ~[2/3  Bm+1/3  B0] 
p  =  coefficient  of  thermal  expansion,  (g/Terg);  NOTE:  BpE  is  equivalent  to 
the  Gruneisen  energy  dependence  of  the  solid  pressure,  i.e.,rpE. 


I 

t 
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3.2. 1.2  Hysteresis  Due  To  Air-Filled  Porosity.  For  a  material  which  contains  air- 
filled  porosity,  i.e.,  less  than  100%  saturation,  the  pores  may  collapse  irreversibly  under 
compression  so  that  the  behavior  of  the  material  is  path  dependent.  Several  additional 
steps  and  a  new  independent  variable  are  required  to  compute  the  pressure.  In  the 
SHEL  routine,  the  new  variable,  pmax  (pm),  is  the  maximum  excess  compression  ever 
attained  by  the  material.  For  a  material  containing  air-filled  void  spaces,  the  initial 
density,  p0,  is  always  less  then  the  reference  density,  i.e.;  po<pr.  By  convention,  the 
excess  compression,  p,  relative  to  p0  is  passed  to  the  SHEL  routine.  Since  the  SHEL 
EOS  algorithms  are  based  on  the  reference  density,  pr,  the  routine  transforms  the  p's  to 
correspond  to  pr. 

A  fundamental  assumption  governs  the  pressure  behavior  of  the  porous  material 
under  loading,  namely,  there  is  a  transformation  in  p  to  an  equivalent  value  along  the 
load  path  of  the  solid  (void-free)  material.  To  calculate  the  pressure,  we  first  find  the 
solid  p  equivalent  to  pm  and  then  obtain  a  modified  p'  for  the  current  compressional 
state  of  the  material  to  use  in  Equation  3.2  as  follows. 

Define: 

'Hz  =  Po  ^  Pr>  Pz  =  Tlz  —  1  —  0>  £z  =  —  Pz  /  T)z  (3-7) 

Any  excess  compression,  p  relative  to  p0>  can  be  transformed  into  an  equivalent 
pr  relative  to  pr,  by 

pr  =  T|zP  +  Pz  (3.8) 

'  As  for  a  nonporous  solid  (Equation  3.4),  we  again  augment  p  and  pm  by  adding 
the  thermal  contribution  to  both; 

p"=p  +  pE 

Pm  =  Pm  +  PE 

An  offset  parameter,  A,  is  computed  next  to  transform  the  compression  into  the 
frame  of  the  reference  density; 

A  =  p2(1  -  ea24m/Ez )  +  p-j(1  - e^m^1  )ea24m/£z  (3  i 

where  a2  and  p1  are  additional  input  parameters.  The  first  term  on  the  right-hand 
side  of  Equation  3.10  provides  a  smooth  convex  curve  starting  at  p=p0  which 


(3.9) 
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asymptotically  approaches  the  reference  P-p  hydrostat;  the  second  term  adds  a  toe  to 
the  initial  loading.  Physically,  the  first  term  represents  the  irreversible  pore  collapse 
and  the  breakdown  of  cementation  bonds  in  the  matrix,  while  the  second  term  is  a 
measure  of  the  initial  elastic  response  of  those  bonds. 

Ps  is  now  computed  as  in  Equation  3.2; 


Ps  =  Bmii-(Bm-BoV'(1-e  “'“"l 

(3.11) 

where 

P  =  TlzP"+A 

(3.12) 

and 

n«=  n  *  /[Tfed-  a2e“2,‘> ) + 

(3.13) 

The  assumption  that  unloading  paths  for  partially  collapsed  porous  media 
parallel  the  solid  load-unload  path  has  been  used  for  many  earth  media,  particularly 
when  load-unload  data  were  not  available  or  inconsistent.  Recommended  load-unload 
paths  for  the  porous  shales  provided  by  WES  did  not  conform  to  this  behavior. 
Unloading  from  peak  pressures  less  than  the  elastic  toe  were  elastic;  as  the  peak 
pressure  increased  above  the  crush  level  so  did  the  initial  unload  moduli.  Thus  the 
unloading  paths  appear  to  stiffen  and  fan  out  as  the  peak  pressure  increases.  This 
effect  was  modeled  by  simply  adjusting  the  effective  maximum  excess  compression  as 
the  material  unloaded,  i.e.,  for  pXi <  pm  <  (xXj 

\ 

pm'=fp  +  (l-f)Pm  (3.14) 

where 

f  =  (M-m  —  M-xi )  ^(M-X2  —  M-xf)  (3.15) 

and  |ixiand  pX2are  input-  Thus,  for  a  porous  solid,  six  additional  parameters  must  be 
specified,  namely; 

p0  =  the  initial  density  of  the  material  (*pr) 
ai  =  relative  slope  of  initial  loading  on  the  elastic  toe  to  that  of  the 
non-porous  solid 
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p ^  =  amount  of  compression  on  loading  toe 

a2  =  fitting  parameter  governing  the  elastic  slope  of  the  load  curve  while 
crushing  out  the  voids 

pxi  =  minimum  compression  to  begin  fanning  the  unload  paths 

p^  =  compression  where  unload  path  parallels  zero  porosity  load-unload  curve. 

3.2. 1.3  Solid-Solid  Phase  Change.  Many  materials  exhibit  a  solid-solid  phase 
change  under  compression.  Based  on  Hugoniot  and  release  data  for  silicates,  it 
appears  that  under  shock  loading  such  a  change  is  hysteretic,  i.e.,  on  unloading  the 
material  initially  remains  in  the  higher  density  state  and  does  not  revert  to  the  original 
state  until  well  down  the  release  path.  Such  behavior  can  be  modeled  by  assuming 
that  the  phase  change  is  a  function  of  both  compression  and  energy. 

To  compute  the  new  solid  PSl  Equation  3.2  or  3.1 1  is  again  modified  by  replacing 
the  constant  Bm  and  the  excess  compression,  p’  (Equation  3.4)  or  p  (Equation  3.12), 
by 


B^  and  p'-A,  or  p-Ap 


(3.16) 


where 


Ap  = 


B'm=  j 


0 

pE  <  (pE)|  (all  state  1) 

9P2 

pE  >  (pE), 

(3.17) 

Bm 

pE  <  (pE).;  (all  statel ) 

(3.18) 

Bm  +  9(Bm2  -  Bm  ) 

pE  >  (pE), 

and 


p2  =  p2/Pr-1. 

1  _  e-«3tPE  ~  (PE)i3 

9=1  +  e-a3[PE-<PE>ll 


(3.19) 


and  Bm2)  (pE)i,  p2  and  0C3  are  material  constants.  Thus,  four  additional  input 
parameters  are  required  for  the  phase  change: 

1  P2  =  the  reference  density  of  the  2nd  phas$ 

'  (pE)i  =  the  energy  per  unit  volume  at  which  the  phase  change  begins 
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«3  =  a  parameter  setting  the  rate  of  conversion  between  the  phases 

Bm2  =  the  maximum  bulk  modulus  (dp/d|i)  of  the  2nd  phase. 

NOTE:  Bm2  is  relative  to  pr  of  the  1st  phase;  the  modulus  relative 

to  the  zero  pressure  density  of  the  2nd  phase,  p2,  is  B2  = 

(p2/pr)  Bm2 

3  2.1.4  Deviatoric  Behavior.  The  complete  description  of  a  material's  stress- 
strain  behavior  requires  the  calculation  of  the  deviatoric  stress  tensor>°  y » in  addition  to 
the  mean  stress  or  pressure. 

The  SHEL  EOS  model  assumes  that  the  stress  deviators  depend  on  the  strain 
deviators,  E'ij,  and  a  yield  function,  Y.  The  stresses  are  calculated  in  two  steps.  First, 
the  change  in  the  deviators  is  assumed  to  be  elastic  and  an  incremental  change  is 
added  to  the  previous  deviators, 

a'jj(n)  =  a'jj(n-1)+2Gde'jj  (3.20) 

where  G  is  the  shear  modulus  of  the  material  and  de'y  is  the  change  in  the  ij-th 
component  of  the  deviatoric  strain  tensor  between  times  n-1  and  n.  Furthermore,  it  is 
assumed  the  shear  modulus  can  be  derived  from  the  current  bulk  modulus  through 
Poisson's  ratio  and  the  standard  elastic  relation 

2G  =  3B(1  -  2v)  /  (1  +  v)  (3.21) 

After  the  stresses  have  been  updated  by  Equation  3.20  a  check  is  made  to  see  if 
the  elastic  changes  are  valid.  The  second  invariant  of  the  deviatoric  stress  tensor,  J'2 
is  compared  to  the  square  of  a  yield  function,  Y  which  depends  only  on  the  mean 
stress,  P.  If  J'2  exceeds  Y2,  the  material  is  deforming  plastically  and  the  deviators  are 
proportionally  reduced  so  that  J'2  is  equal  to  Y2.  For  the  rock  media  of  interest,  the 
yield  function  is 

-P/PV 

Y  -  Yvm  -  (Yvm  -  Y0)e 


(3.22) 


The  four  input  parameters  required  to  calculate  the  deviatoric  stress  tensor  are; 
v,  Poisson's  ratio 

Yvm,  the  limiting  yield  strength  (kb) 

Y0,  the  yield  strength  at  zero  pressure  (kb) 

Py  the  pressure  at  which  the  yield  function  —(2/3  Y^  +  1/3  Y0)-  (kb) 


Table  3-5.  Equation-of-state  parameters  for  SHEL  (EOS  1). 


Value 

Units 

Name  CK  Location 

Comments  /  Definition 

Po 

g/cc 

RHOZ 

Initial  density 

Pr 

g/cc 

RHOREF 

Reference  density 

Es 

Terg/g 

ES 

Min  energy  for  Pv 

Bo 

Mb 

BO 

Initial  bulk  modulus 

Bm 

Mb 

BM 

Maximum  bulk  modulus 

* 

P 

— 

XMU 

Rate  of  change  of  B 

PI 

— 

End  of  elastic  toe 

ai 

— 

Rate  of  increase  of  B  on  elastic  toe 

a2 

— 

Rate  of  increase  of  B  during  pore  collapse 

Pxl 

— 

Start  of  unload  path  fan 

Px2 

— 

End  of  unload  path  fan 

P2 

g/cc 

Ref  density  of  higher  density  phase 

(pE)l 

Te/cc 

Min  E N  for  solid  in  2nd  phase 

Bm2 

Mb 

Max  B  of  2nd  phase  solid 

— 

Rate  for  solid-solid  phase  change 

Vv 

% 

Porosity  (total  void  volume) 

P  - 

g/Terg 

CTE 

Thermal  expansion  coef. 

V 

— 

Poisson's  ratio 

Yvm 

Mb 

Von  Mises  yield  limit 

Yo 

Mb 

Cohesion 

Pv 

Mb 

Coef.  in  Rate  of  increase  in  Y  with  P 

20 


3.2.2  VIC. 


Table  3-6.  Equation-of-state  parameters  for  VIC  (EOS  2). 


VALUE 

J/M 

UNITS 

1 

A 

1 

— 

2 

B 

1 

- 

3 

EO 

1 

Te/g 

4 

ESS 

1 

Te/g 

5 

ESP 

1 

Te/g 

6 

ALPHA 

1 

-- 

7 

BETA 

1 

— 

8 

AA 

1 

Mbar 

9 

BB 

1 

Mbar 

10 

GM 

1 

Mbar 

11 

POIS 

1 

— 

12 

GOK 

1 

-- 

13 

DELE 

1 

14-16 

17 

AP 

1 

Mbar 

18 

BP 

1 

Mbar 

19 

EPM 

1 

Te/g 

20 

EPP 

1 

Te/g 

21 

DEP 

1 

— 

22-27 

28 

BETAH 

1 

29 

BETAS 

1 

-- 

30 

EPYLD 

1 

— 

31 

YFLOWR 

1 

32 

YLDMXF 

1 

Mbar 

33 

YLDMXC 

1 

Mbar 

»34 

COHC 

J 

Mbar 

35 

DYDPC 

1 

— 

36 

YLDMXE 

1 

Mbar 

37 

COHE 

1 

Mbar 

38 

DYDPE 

1 

— 

39 

COHF 

1 

Mbar 

40 

DYDPF 

1 

— 

41-48 

49 

VREF 

cm3/g 

50 

NEOS 

2- 

DEFINITION 
y-1  of  vapor  phase 
Coef  of  2nd  vapor  term 
Parameter  in  vapor  term 
Energy  of  incipient  vaporization 
Energy  of  complete  vaporization 
Exponential  decay  coef 
Exponential  decay  coef 
Bulk  modulus  at  zero  P 
Coefficient  of  jx2  term 
Max  Shear  modulus 
Poisson's  Ratio,  y 
G/K  [=1 .5(1  -2v/(1+v) 

1/(ESP-ESS) 

AA  after  phase  change 

BB  after  phase  change 

Min.  E  for  solid-solid  phase  change 

Max.  E  for  solid-solid  phase  change 

_.  of  phase  change 

Rate  of  hardening  of  Y 

Rate  of  softening  of  Y 

Plastic  |x  where  hardening  ends 

Flow  rule  sentinel  (0=non-assoc.,  1  =  assoc.) 

Von  Mises  limit  if  fractured  (def=YLDMXC) 

Von  Mises  yield  limit,  (J3’>0) 

Cohesion  of  competent  material  (J3'>0) 
dY/dP  in  compression  (J3'>0) 

Ymax  in  extension  (J3' 

Cohesion  in  extension  (J3' 
dY/dP  in  extension  (J3‘ 

Cohesion  of  3-cracked  material  (def=COHC) 
dY/dP  after  3-cracks  (def=DYDPC) 

1/RHOREF 
EOS  #  for  VIC 


I 
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3.2.3  CIST. 


Subroutine  CIST  is  a  relatively  simple  EOS  developed  at  the  AFWL5  to  model 
the  results  of  Cylindrical  In-Situ  Tests  (CIST)  field  data.  It  has  been  used  for  both  low 
and  high  stress  calculations,  such  as  those  in  Reference  12.  The  current  CRALE 
version  uses  the  AFWL  model  for  the  pressure-density  relation,  but  calculates  the 
energy  dependence  in  the  same  way  as  the  SHEL  EOS.  Low  pressure  loading  and 
unloading  occur  along  straight  line  segments  determined  by  the  sets  of  input  wave 
speeds,  Cl,  C2,  C3,  pressure  limits,  PI,  P2,  P3,  and  Poisson's  ratios  PRL  and  PRU. 
For  pressures  greater  than  P3,  the  routine  uses  a  quadratic  similar  to  the  SHEL  EOS, 
i.e., 


P  =  P3+B1(|i-p3)  +  Bsq(p-p3)2  (323) 

At  high  energy  densities,  an  additional  term  is  added  to  the  pressure.  This  term 
is  calculated  in  the  same  way  as  in  subroutine  SHEL  (Equations  3.1  through  3.4). 


I 

I 

5  AFWL  ref 
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Table  3-7.  Equation-of-state  parameters  for  CIST  (EOS  3). 


£K 

VALUE 

I/M 

UNITS 

1 

Cl 

Y 

cm/ms 

2 

C2 

Y 

cm/ms 

3 

C3 

Y 

cm/ms 

4 

PI 

Y 

Mbar 

5 

P2 

Y 

Mbar 

6 

P3 

Y 

Mbar 

7 

BSQ 

Y 

Mbar 

8 

ALFA 

Y 

9 

EMU3 

Y 

--  - 

10 

EMU2 

Y 

— 

11 

EMU1 

Y 

— 

12 

EMUMX 

Y 

— 

13 

B1F 

M 

14 

B1 

Y 

Mbar 

15 

B2 

Y 

Mbar 

16 

B3 

Y 

Mbar 

17 

CTE 

Y 

g/Te 

18 

EV 

Y 

Te/g 

19 

FPRL 

Y 

— 

20 

FPRU 

Y 

— 

21-22 

23 

BETA 

Y 

24 

25 

PRL 

Y 

26 

PRU 

Y 

— 

27 

GMAX 

Y 

Mbar 

28 

BETAH 

1 

— 

'  29 

BETAS 

1 

— 

x  30 

EPYLD 

— 

31 

YFLOWR 

'l 

32 

YLDMXF 

1 

Mbar 

33 

YLDMXC 

1 

Mbar 

34 

COHC 

1 

Mbar 

35 

DYDPC 

1 

— 

36 

YLDMXE 

1 

Mbar 

37 

COHE 

1 

Mbar 

38 

DYDPE 

1 

— 

39 

COHF 

1 

Mbar 

40 

DYDPF 

1 

— 

41-48 

49 

VREF 

Y 

cm3/g 

50 

NEOS 

Y 

3 

DEFINITION 

Unload  and  initial  load  wave  speed 
2nd  loading  wave  speed 
3rd  loading  wave  speed 
P  where  C  changes  from  Cl  to  C2 
P  where  C  changes  from  C2  to  C3 
P  where  load/unload  curves  merge 
Coef.  of  p  in  high  pressure  fit 
Controls  slope  of  unloading  heel 
p  where  load/unload  curves  merge 
pat  P2 
pat  PI 

p  at  P=0  along  max  unload  curve 

Ratio  of  initial  K  to  B1  (def=1 .) 

dP/dp  corresponding  to  Cl 

dP/dp  corresponding  to  C2 

dP/dp  corresponding  to  C3 

P,  the  Coefficient  of  thermal  expansion 

Min.  vaporization  energy  (>  ES) 

f|_  =  1 .5(1  -2v|_)/(1  +vl)  =  G/K  [loading] 

fu  =  1 .5(1  -2vu/(1  +vu)  =  G/K  [unloading] 

Controls  break  on  unloading  for  heel 

Loading  Poisson's  ratio,  vl 

Unloading  Poisson's  ratio,  vu 

Maximum  shear  modulus 

Rate  of  hardening  of  Y 

Rate  of  softening  of  Y 

Plastic  p  where  hardening  ends 

Flow  rule  sentinel  (0=non-assoc.,  1=  assoc.) 

Von  Mises  limit  if  fractured  (def=YLDMXC) 

Von  Mises  yield  limit,  (J3'>0) 

Cohesion  of  competent  material  (J3'>0) 
dY/dP  in  compression  (J3’>0) 

Ymax  in  extension  (J3'<0) 

Cohesion  in  extension  (J3'<0) 
dY/dP  in  extension  (J3'<0) 

Cohesion  of  3-cracked  material  (def=COHC) 
dY/dP  after  3-cracks  (def=DYDPC) 

1/RHOREF 

the  material  uses  the  CIST  EOS 


I 

I 
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3.2.4  HE. 


A  model  for  High  Explosives,  air,  simple  gases  or  mixtures  of  HE  and  air. 

HE  is  a  combination  of  the  Lawrence  Livermore  National  Laboratory  (LLNL)  JWL6 
model  for  high  explosives  and  an  AFWL7  variable-gamma  or  constant-gamma, 
perfect-gas  model  of  air.  For  HE,  the  JWL  model  calculates  pressure  as  the  sum  of 
three  terms,  a  gamma-law  gas  plus  two  terms  which  decay  exponentially  with  increasing 
volume;  namely; 

P  =  A  *  (1  -  to  /  R1v)e-R1v  +  B  *  (1  -  ©  /  R2v)e_R2v  +  copE  (3  24) 

Values  of  the  JWL  parameters  for  several  typical  HEs  are  given  in  Table  3-8  and 
defined  in  Table  3-9.  The  input  data  for  the  JWL  model  are  stored  in  CK  1-8.  Note  the 
input  units  for  E0,  the  initial  explosive  energy  [CK(7)]  are  Te/cc  rather  than  the  code 
units  of  Te/gm.  This  choice  was  made  to  retain  consistency  with  general  usage.  If  CK1 
is  0.,  the  material  is  assumed  to  be  air  or  a  perfect  gas;  Equation.  3.24  is  skipped  and 
the  air  or  perfect  gas  algorithm  is  used. 

The  release  adiabat  for  TNT,  shown  in  Figure  3-1 ,  illustrates  the  effect  on  the 
total  pressure  of  each  term  on  the  right-hand  side  (rhs)  of  Equation  3.24.  Both  the  first 
(A-term)  and  second  (B-term)  parts  of  the  rhs  of  Equation  3.24  decay  very  rapidly  as  the 
HE  products  expand.  As  shown  in  Figure  3-1,  both  terms  become  insignificant  relative 
to  copE  at  volumes  greater  than  about  5  cc/g,  well  above  the  densities  usually  attained  in 
'  air  or  similar  gases.  Hence  at  large  volumes  the  JWL  model  reduces  to  a  simple, 
constant-y  law  gas,  where  co=y-1 .  Because  many  problems  of  interest  involve  HE 
explosions  in  air,  modeled  as  a  variable-y  gas,  the  JWL  model  has  been  modified  so 
that  at  large  volumes  the  effective  y  transions  from  the  JWL  value  to  that  calculated  by 
the  AFWL  algorithms. 

The  transition  from  the  JWL  gamma  to  that  of  air  begins  at  the  density  where  the 
B-term  is  less  than  PCUT,  currently  set  to  10‘3  bars,  i.e.; 

Pb  =  R2po/loge[PCUT/B]  (3.25) 


I 

I  _ 

®  LLL  HE  reference 
7  AFWL  ref  for  air  EOS 
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and  and  is  equal  to  the  air  y  for  densities  less  than  or  equal  to  .01  g/cc,  about  a  factor 
of  ten  greater  than  ambient  air.  A  linear  interpolation  in  density  is  used  to  go  smoothly 
from  one  oo  to  the  other.  By  allowing  co  to  vary  as  a  function  of  the  gas  density,  a  single 
EOS  may  be  used  for  pure  HE,  pure  air,  or  a  mixture  of  the  two.  This  is  convenient 
when  calculating  an  HE  explosion  in  air  since  the  boundary  between  the  two  materials 
may  then  be  treated  as  non-Lagrangian. 

If  AGAM  is  set  to  zero,  the  routine  assumes  the  material  is  pure  explosive;  if 
AGAM  is  greater  than  zero,  the  material  will  be  calculated  either  as  pure  air  [HEGAM  = 
0]  or  a  mixture  of  explosive  and  air  ( HEGAM  >  0).  If  AGAM  is  not  zero,  EOS  HE  will 
use  the  AFWL  variable  gamma  model  for  the  gas.  To  model  a  constant  gamma-law 
gas,  set  AGAM  to  1 HEGAM  to  0.,  and  GAMPG  to  the  desired  value  of  y-1 . 


I 

t 
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PRESSURE  (Mb) 


Figure  3-1 .  Release  Path  from  the  CJ  point  of  TNT  showing  the 
contributions  of  the  three  terms  in  the  JWL  EOS. 
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Table  3-8.  JWL  parameters  for  some  common  high  explosives. 


HE 

Po 

PCJ 

Eo 

D 

(0 

A 

R1 

B 

R2 

type 

g/CC 

kb 

Te/cc 

cm/»s 

(ri) 

Mb  ' 

Mb 

TNT 

1.63 

210 

.06 

.693 

.35 

3.738 

4.15 

.03747 

0.90 

NM 

1.128 

125 

.051 

.628 

.30 

2.092 

4.4 

.05689 

1.2 

NM1 

1.139 

125 

.050 

.6187 

.27 

1.922 

4.3 

.0631 

1.3 

ANFO 

0.82 

55 

.0315 

.455 

.31 

0.4716 

3.7 

.009536 

0.95 

PETN 

1.77 

335 

.101 

.83 

.25 

6.17 

4.4 

.16926 

1.2 

9404 

1.84 

370 

.102 

.88 

.25 

8.545 

4.6 

.20493 

1.35 

Table  3-9.  Equation-of-state  parameters  for  HE  (EOS  4). 


muz 

m 

UNITS 

DEFINITION 

1 

HEGAM 

i 

to  =  y-1  for  High  Explosive 

2 

A 

i 

Mb 

parameters  of  the 

3 

R1 

i 

J-W-L  (LLL)  equation 

4 

B 

i 

Mb 

of  state,  Equation  3.24 

5 

R2 

i 

(set  to=  0  for  EOS  of  air) 

6 

DETV 

i 

cm/ps 

Detonation  velocity 

7 

EPOT 

i 

Te/cm3 

Potential  energy  released  on  detonation 
(divide  by  p0  to  get  specific  energy,  EZ) 

8 

CJP 

i 

Mbar 

Chapman-Jouguet  Pressure 

9 

AGAM 

i 

=  0  for  HE  only;  =  1  for  air  or  mix  of  air-HE 

10 

GAMPG 

i 

if  a)  =  0.,  value  of  y-1  for  perfect  gas 

11-16 

not  used 

17 

■A1 

M 

g/cm3 

Coef.  of  p  in  A  term,  Equation  3.24 

18 

A2 

M 

cm3/a 

HEGAM/A1 

19 

CUTA 

M 

g/cm3 

-A1/LOG(10'9/A) 

20 

B1 

M 

g/cm3 

Coef.  of  p  in  B  term,  Equation  3.24 

21 

B2 

M 

cm3/a 

HEGAM/B1 

22 

CUTB 

M 

g/cm3 

-B1/LOG(10'9/B) 

23 

24-30 

CUTC 

M 

g/cm3 

.01 ,  lower  density  limit  for  HE  /  air  mixture 
not  used 

31 

CJMU 

M 

excess  compression,  p,  at  CJ  point 

32 

33-48 

CJGAM 

M 

cm3/gm 

effective  7-1  at  CJ  point 
not  used 

49 

VREF 

M 

1/RHOREF 

50 

NEOS 

M 

4. 

-  the  material  uses  the  HE  EOS 

I 
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3.2.5  MCST. 


A  multi-segmented  "stick"  model  based  on  the  CIST  model  has  been 
constructed,  primatily  to  fit  the  behavior  of  the  Socorro  plaster  sand  used  in  numerous 
tests  at  the  DNA  PHETS  site  at  White  Sands.  The  model  is  quite  flexible  and  is  also 
used  to  model  concrete.  For  details  contact  Jim  Rocco  at  the  TRT  in  Albuquerque. 

3.2.6  SCUB. 

This  module  is  currently  inactive.  It  has  served  in  the  past  as  the  location  for  the 
table  lookup  routines  provided  by  SCUBED  for  material  models  used  in  their  early  time 
cratering  calculations  in  several  DNA  programs,  e.g.,  Benchmark  Cratering  and  MISTY 
ECHO.  Since  SCUBEDs  EOS  models  are  also  dynamic,  it  is  reccommended  that  they 
be  contacted  for  their  current  model  if  it  is  to  be  used  in  CRALE  or  EQS  calculations. 

3.2.7  H20. 

H20  contains  the  SAIC  1983  model  for  water.  This  model  was  incorporated  to 
compare  with  others  and  is  not  recommended  for  use  without  further  study. 

3.2.8  UEOS. 

UEOS  is  an  empty  routine  included  to  provide  the  hooks  by  which  users  may 
easily  incorporate  their  own  EOS  model  into  CRALE. 


I 

I 
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3.2.9  CAP. 


The  CAP  EOS  currently  in  use  in  CRALE  was  developed  by  Dr.  George  Baladi 
while  at  the  Waterways  Experiment  Station  (WES).  It  is  a  strain  dependent  model  in 
that  the  complete  stress  tensor  is  updated  each  cycle  based  on  the  current  stress  state 
and  the  incremental  change  in  strain.  The  model  also  calculates  the  shear  modulus,  G, 
at  the  average  strains  rather  than  using  an  initial  or  final  value  as  is  required  when  G  is 
a  function  of  P.  A  problem  with  previous  CAP  models  was  their  sensitivity  to  step  size. 
Baladi's  model  reduced  this  condition  by  using  the  average  moduli  on  each  cycle. 

Table  3-10.  Equation-of-state  parameters  for  CAP  (EOS  9). 


£K  VALUE  UM  IMIS  DEFINITION 


1 

AKEI 

Y 

MPa 

2 

AK1 

Y 

ft 

3 

AK2 

Y 

N 

4 

AK3 

Y 

ft 

5 

AK4 

Y 

ft 

6 

AK5 

Y 

It 

7 

AK6 

Y 

N 

8 

AK7 

Y 

N 

9 

AKIM 

Y 

10 

AKIM 

Y 

It 

11 

AK2M 

Y 

H 

12 

AA 

Y 

13 

AB 

Y 

14 

AB1 

Y 

It 

15 

AC 

Y 

17 

AGEI 

Y 

18 

AG1 

Y 

tt 

19 

AG2  , 

Y. 

II 

20 

AG3 

Y 

H 

21 

AG4 

Y 

It 

22 

CRO 

Y 

23 

CR1 

Y 

It 

24 

CR2 

Y 

It 

25 

AW 

Y 

26 

AD 

Y 

27 

ADI 

Y 

It 

28 

AD2 

Y 

tt 

29 

AD3 

Y 

tt 

30 

AD4 

Y 

tl 

33 

SJMX 

Y 

34 

CKMX 

Y 

35 

CKAA 

Y 

36 

LTYPE 

Y 

41 

ASTART 

42 

XSTART 

43 

FCUT 

44 

TCUT 

Used  to  calculate  K 

It  ft  M 

ft  M  M 

H  ft  H 

N  ft  N 

ft  M  H 

It  tt  It 

It  ft  It 

Used  for  pore  pressure  K 

tt  h  tt  it 

tt  N  N  M 

Maximum  Yield  stress 
Coef  in  exp  change  of  Y 

tt  It  M  It  It  tt 

Cohesion  =  AA-CC 
Used  to  calc  G  (Shear  Mod) 

tt  tt  tt  it  tt  it 

H  It  It  H  II  II 

M  It  ft  ft  tt  M 

H  H  It  M  It  It 

Cap  parameter 

It  H 

H  II 

Maximum  hysteretic  strain 
Used  to  calc  hysteretic  strain 

tt  it  it  tt  it 

it  it  it  tt  n 

M  ft  It  It  It 

It  tt  H  It  M 

Pressure  switch  to  hi-P  EOS 
Maximum  bulk  modulus 
exp  rate  of  increase  to  CKmx 
=1 ,  cap  retracts;  =2,  fixed 
Initial  position  of  L 
Initial  Position  of  X 
Max  Tensile  pressure 
Tensile  cutoff 
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3.2.10  MIX. 


Table  3-1 1 .  Equation-of-state  parameters  for  MIX  (EOS  1 0).  (Table  to  be 
included  in  next  issue  of  manual) 

3.2.11  Others. 

Over  the  years  of  using  AFTON  and  CRALE  a  number  of  EOS  models  have 
been  incorporated  in  various  PL's.  Several  may  still  be  included  on  the  current  PL,  but 
their  use  would  require  some  preliminary  testing  to  insure  that  they  still  are  operational. 
They  are  included  here  for  completeness. 


3.2.1 1.1  TILL.  TILL  is  a  modified  version  of  the  TILLOTSON  EOS8.  This  EOS 
was  derived  originally  for  very  high  pressures  in  metals.  For  densities  greater  than  the 
initial  RHOZ,  the  pressure  is  calculated  as 


P  =  A|i  +  Bp2  + 


a-i- 


b/(E/E0n2  +  l)]pE 


(3.26) 


If  a  Poisson's  ratio  (POIS)  and  nonzero  value  of  Gm,  and  the  maximum 
allowable  shear  modulus  are  specified,  TILL  will  return  the  elastic  moduli  to  STRAIN  so 
deviatoric  stresses  may  be  computed.  The  TILL  EOS  will  not,  however,  treat  hysteretic 
or  plastic  materials.  Although  it  has  been  used  for  rocks,  the  low  stress  regime  (  CIST). 


3.2.1 1 .2  GRAY.  The  GRAY9  EOS  is  a  version  of  the  routine  developed  by  LLNL 
'  for  metals.  It  does  not  have  deviatoric  components  and  may  give  erroneous  results  at 
'very  high  pressures. 


I _ 

®  J.  Tillotson,  "Metallic  Equations  of  State  for  Hypervelocity  Impact",  GA-3216,  General  Atomic,  July  1962. 
9  E.B.  Royce,  "GRAY,  A  Three-Phase  Equation  of  State  for  Metals",  UCRL-51121,  LLL,  Sept.  1971. 
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SECTION  4 


EQS  -  THE  EQUATION  OF  STATE  DRIVER 

4.1  OVERVIEW. 

EQS  is  a  stand-alone  program  designed  to  exercise  the  EOS  models  described 
in  the  previous  sections  that  are  used  in  the  1-  and  2-D  CRALE  codes.  The  program 
has  a  number  of  options  for  both  specific  and  arbitrary  stress-strain  paths  and  generates 
plots  for  easy  comparisons  with  various  types  of  material  properties  data. 

An  important  feature  of  EQS  is  its  conformity  with  the  CRALE  codes.  Input  data 
for  the  various  EOS  models  are  read  using  EQST,  the  same  routine  called  by  CRALE  1 
and  CRALE2.  The  major  subroutines  in  EQS  generate  the  desired  stress-strain  paths 
by  first  setting  up  an  axisymmetric  unit  cell  in  the  2D  code,  prescribing  the  appropriate 
motions  of  that  cell's  corners,  and  then  calling  the  CRALE2  subroutine  STRAIN  to 
calculate  the  incremental  strains  and  resulting  stresses.  Thus,  once  a  model  has  been 
checked  out  with  the  driver,  the  same  behavior  should  be  followed  in  the  ID  and  2D 
codes.  (However,  note  that  since  EQS  only  exercises  a  single  j,k  cell  in  a  pseudo-2D 
CRALE  grid,  subtle  errors  in  setting  local  variables  could  still  occur  in  the  FD  codes.) 

The  current  EQS  program  stores  up  to  200  complete  sets  of  stress-strain  data  (a, 
e,  p,  E,  moduli,  speeds,  etc.)  for  each  Hugoniot,  release  adiabat  or  strain  path  generated 
'  by  HUG,  REL,  PATH  or  FLIP  and  experimental  data  (a,  p,  Up,  Us,  E)  for  each  call  to 
EXP  or  EQN.  Each  data  set  is  identified  by  a  line  number  [LINES]  and  is  referred  to  as 
such  when  needed  by  subsequent  calls  to  other  routines.  After  the  data  sets  have  been 
produced,  selecting  the  PLOT  option  provides  the  flexibility  to  overlay  data  for  different 
EOS’s  and/or  the  experimental  Hugoniot  data  or  fit  on  a  single  plot.  The  data  may  also 
be  written  to  an  output  file  for  further  analysis  by  other  programs. 


I 
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4.2  CODE  MODULES. 


4.2.1  EQS,  the  Main  Program. 

The  main  program,  EQS,  serves  as  the  interface  between  the  user  and  the 
various  routines  that  generate  the  data  sets  and  plots.  The  basic  flow  of  the  EQS 
program  is  illustrated  by  the  schematic  flowchart  in  Figure  4-1 ;  the  input  file  required  to 
run  EQS  is  summarized  in  Table  4-1 .  Details  of  the  input  as  used  by  each  subroutine 
are  described  in  the  subsections  detailing  that  routine.  EQS  is  quite  fast,  typically 
completing  a  run  in  under  30  seconds  on  a  486  PC.  Hence,  the  remainder  of  this 
section  assumes  the  operations  occur  on  a  PC.  The  code  will  run  in  batch  mode  on  the 
DNA  workstations  and  CRAY  computers  with  only  minor  differences. 


EQS  begins  by  calls  to  several  subroutines  that  initialize  the  run.  The  first  call  to 
OPNFIL,  opens  the  appropriate  I/O  and  data  files.  Currently  the  various  files  that  may 
be  required  during  execution  are; 
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Unit  No. 

2 

3 

5 

6 
10 
11 

12  ■ 

14 


PC  File  Name 

’name'.ZTA 

EOSS.DAT 

’name’.INP 

'name'. OUT 

HUGROCK.DAT 

HUGMET.DAT 

RELAD.TXT 

USUP.DAT 


Purpose 

contains  the  plot  output 

A  data  base  of  specific  material  EOS  input 

The  input  file  read  by  EQS 

The  standard  output  file 

Experimental  Shock  data  for  rocks,  grouts,  H20 

Experimental  Shock  data  for  metals,  plastics,  etc. 

Release  Adiabat  data  from  SNL  Hugoniot  points 

Data  base  of  fits  to  US-UD  on  the  Hugoniot. 


where  'name'  implies  a  file  name  supplied  by  the  user. 


EQS  obtains  its  operating  instructions  from  an  ASCII  input  stream  constructed  by 
the  user  as  a  name.inp  file.  Line  1  of  the  input  (see  Table  4.1 )  includes  a  user  ID 
(PLTNAM),a  switch  ( /CA/V )  to  set  the  I/O  units  (default  =  Mb-g-ps)  and  contols  for 
including  a  logo  on  the  plot  output  (NLOGO,XLPLT,  YLPLT,SIZEL).  The  code 


^  Although  the  ICNV  parameter  may  be  used  to  set  the  I/O  units,  some  output,  particularly  when  printed  using  the 
'  F-type  format,  remain  in  code  units  (Mb-g-ps).  Wave  speeds  and  moduli  in  HUG,  for  example,  are  printed  in 
code  units,  i.e.,  cm/ps  and  Mb. 
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continues  with  a  call  to  INITIAL,  a  general  CRALE  routine  that  initializes  parameters, 
such  as  k  and  the  I/O  units  to  be  used  in  the  run.  Although  eight  sets  of  units  are 
available,  if  ICNV  is  set  to  9  or  10,  INITIAL  will  read  in  two  lines  containing  the  desired 
conversion  factors  from  code  to  I/O  units  and  corresponding  names.  While  the  code 
uses  a  consistent  set  of  units  in  its  calculations,  these  I/O  units  may  be  completely 
arbitrary.  The  CRALE  utilities  library  contains  a  routine  to  draw  the  TRT  logo.  Other 
logos  can  be  added  to  the  library  and  then  placed  on  EQS  plots  by  setting  NLOGO 
greater  than  zero.  XLOGO,  YLOGO  and  SIZEL  can  be  set  to  position  the  logo  and 
adjust  its  size  from  the  defaults. 


line  1 

PLTNAM 

ICNV 

NLOGO 

XLOGO 

YLOGO 

SIZEL 

lins.2 

line  3 
NEOS 

TITLE  ' 

line  4 
ROUT 
NP 

SETUP 


Table  4-1 .  Summary  of  the  input  to  EQS. 

FORMAT  (A4,  212,  12X,  3F5.0) 

4  Character  ID  for  plots 

Specifies  I/O  units;  default  is  Mb-gm-ps 

0=  omit  logo  from  plots;  1  =plot  TRT  logo;  2=? 

X  position  of  logo  on  page,  {def=  left  edge  of  page} 

Y  position  of  logo  on  page,  {def=  2'  from  top  of  page} 

Size  of  logo,  {def=.5"} 

EOS  data  input;  see  Section  2.5 

FORMAT  (A8,  2X,  18A4) 

Name  of  EOS  to  be  exercised;  blank  terminates  the  run. 
’DATA'  may  be  used  for  NEOS  for  calls  to  PLOT,  EQN  or  EXP. 
Title  for  succeeding  print  and  plot  output 

FORMAT  ( A4,  812,  12F5.0)  see  Table  4.2 
Name  of  EQS  subroutine  to  be  exercised 
8  integer  parameters  for  the  various  subroutines 
Up  to  12  floating  point  parameters  for  the  various  subroutines 


line  5  Additional  data  lines  read  in  EQN  or  PATH 


At  this  point  EQS  is  ready  to  process  the  user  input.  The  program  continues 
with  calls  to  two  entry  points  in  the  EQST  subroutine,  ESINIT  which  reads  in  the 
parameters  of  the  EOS  models  required  for  the  run  (lines  2)  and  ESPRNT  which 
echoes  them  to  the  name. out  file.  ESINIT  will  read  data  for  up  to  15  different  materials, 
each  uniquely  identified  by  name  [IFO$(\)].  (At  least  material  model  must  be  read.) 
Yhe  EQST  input  was  discussed  in  Section  3.1  and  is  summarized  in  Tables  3-1  and 
3-2. 
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After  the  material  models  are  read,  EQS  can  finally  cycle  through  the  sets  of 
user  requests  for  information.  Each  set  of  input  data  begins  with  line  3  containing  the 
name  of  the  material  ( NEOS )  to  be  acted  on  and  an  optional  title  ( TITLE)  for  print  and 
plot  output.  The  code  checks  to  be  sure  the  material  requested  was  one  of  the  EOS 
models  read  in.  (The  word  DATA  will  also  be  accepted  as  a  valid  NEOS  to  process  the 
experimental  data  base.)  If  the  requested  NEOS  is  not  available,  the  code  prints  an 
error  message  and  aborts.  For  a  valid  NEOS,  EQS  reads  the  input  on  line  4  containing 
a  router  {ROUT),  and  as  many  as  8  integer  {NP)  and  12  floating  point  {SETUP) 
parameters. 

ROUT  determines  which  subroutines  to  activate  and  NP  and  SETUP  activate  the 
various  options  and  set  the  parameters  required  by  that  routine  as  summarized  in 
Table  4-2.  Several  of  the  subroutines  (e.g.,  EQN  ,  PATH)  may  ask  for  additional 
input,  line  5,  which  is  read  in  next.  The  formats  for  these  data  are  included  in  the 
subsections  describing  those  routines.  After  completing  the  desired  operations,  the 
code  loops  back  to  read  another  line  4  and  performs  the  next  request.  A  blank  line  4, 
i.e.,  ROUT-  ',  causes  the  code  to  loop  back  to  read  another  line  3  with  the  next 
material  model  to  process.  The  end  of  the  run  is  signaled  by  a  blank  line  3,  i.e., 

NEOS-  the  code  closes  all  open  files  and  terminates.  The  current  choices  for 
ROUT  are; 


ROUT  = 

pouting 

Sect 

Function 

'  EQN 

\ 

EQN 

4?2.2 

generate  the  Hugoniot  data  set  using  US-UD  or  P-UD  relation 
from  least-squares  fit  to  EXP  data  set  or  database 

EXP 

exphug’ 

4.2.3 

Retrieve  experimental  Hugoniot  and  release  adiabats  from 
the  DNA  shock  database 

FLIP 

FLIP 

4.2.4 

Flip  and  translate  Hugoniots  to  get  shock  interactions  and 
BCFs2 

HUG 

HUG 

4.2.5 

Generate  principal  Hugoniot  for  selected  EOS  model 

PATH 

PATH 

4.2.7 

Generate  Cjj  vs  ejj  for  arbitrary  stress-  or  strain-driven  paths 

PLOT 

PLOTN 

4.2.8 

Plot  previously  data  sets  in  several  standard  forms 

REL 

blank 

REL 

4.2.6 

Generate  release  adiabats  from  points  along  the  last 
Hugoniot  or  exercise  model  for  specific  p.,E  pairs. 

Close  files  and  terminate  run 

A  more  complete  description  of  the  NP  and  SETUP  input  parameters  is  included 
in  the  discussion  of  the  individual  routines,  Section  4.2.2  through  4.2.8. 


2  Borehole  Correction  Factors 


Table  4-2.  Parameters  on  line  4  of  EQS  input. 
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4.2.2  EQN. 


For  ROUT=EQN,  subroutine  EQN  is  selected  to  produce  a  Hugoniot  derived 
either  by  fitting  a  previous  data  set  either  1)  generated  by  EXPHUG  [NP(2)=<1]  or 
2)from  input  of  the  coefficients  of  a  quadratic  fit  to  Us-Up  [NP(2)=2]  or  P-Up  [NP(2)=3] 
curves  or  3)  from  a  fit  stored  in  the  Us-Up  DNA  data  base  on  Unit  14  [NP(2)=4] .  The 
Input  to  EQN,  summarized  in  Table  4.2,  is  as  follows; 


HE 

1  LINE  User  Reference  Number  of  data  set  for  plotting,  etc 

2  IEQ  =  0, 1 ;  obtain  least  squares  fit  (LSF)  of  Us-Up  data  for  line  NP(3) 

2;  generate  Hugoniot  from  LSF  to  Us-Up  input 

3;  "  "  . .  P-Up  input 

4;  read  Us-Up  fit  from  USUP.dat  data  base  (Unit  14) 

3  LIN  identifier  of  previous  data  set  to  be  fit,  needed  if  HEI21  ^  *1 


8 

NPP  if  >  0,  line  is  segmented,  read  additional  data  (line  5)  see  below 

forH£I21<1 

get  least  squares  fit  for  LIN;  generate  Hugoniot  data  set 

SETUP 

1 

RHOO 

Initial  density,  only  J=1  is  needed 

2 

UMN 

Minimum  Up  on  segment 

3 

UMX 

Maximum  Up  on  segment 

if  NPP  >  0;  read  in  one  additional  line  (IX,  F9.0,  7F10.0)  with  UMX(i),i=1,8 

> 

for  NPr21  =  2. 

3  read  in  Us-Up  or  P-Up  relation  and  generate  data  set 

SEIUE 

i  =  1 

RHOO 

Initial  density,  only  input  for  J=1  is  needed 

2-4 

Ci 

where  <I>  =  Cq  +  C-|Up  +  C2Up2  and 

<X>  =  Us  if  NP[1]=2;  P  if  NP[1]=3. 

5 

UPMIN 

Minimum  Up  on  segment  {0.0001  cm/ps} 

6 

PMAX 

Maximum  UD  or  P  on  segment  if  NP(1)=2  or  3  {50  Mb} 

if  NPP  >  0;  read  additional  lines  (11,  F9.0,  5F10.0)  containing 

NPP, Cl, C2.C3, UPMIN, PMAX  until  NPP=0. 

for  NPf21  =  4 

read  in  Us-Up  from  data  base  on  Tape  14 

SETUP 

1 

ENUM 

Material  number  in  data  base  (See  Table  5-3.) 

'  2 

IAUT 

Author  ID 

•  3 

RHOO 

Initial  density 
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NP[21<1 :  EQN  finds  the  data  set  corresponding  to  LIN=NP[3]  and  does  a  least 
squares  fit  to  find  the  linear  Us-Up  relation, 

Us  =  C0  +  slip  (4.1) 

and  the  standard  deviation,  a,  of  the  fit.  The  Us-Up  data  may  be  fit  in  up  to  8  intervals. 

If  NPP  *  0,  EQN  reads  a  line  of  data  (8F10.0)  containing  the  UMX  of  each  interval 
[The  UMX  of  each  interval  is  the  default  UMN  for  the  next,  the  initial  UMN=SETUP(2).] 
The  routine  generates  and  saves  the  full  set  of  Hugoniot  parameters,  including  dP/dp 
(defined  as  [dP/dUp]/[Dji/dUp]).  The  output  includes  a  comparison  of  the  shock  P  and 
Us  calculated  by  the  fit  with  the  those  in  the  data  set  used  for  the  fit.  EQN  then 
removes  all  data  that  are  outside  the  ±3ct  bounds  of  the  fit  and  refits  the  remaining 
points.  This  fit  is  not  saved  but  a  printout  comparing  it  to  the  data  is  also  generated. 

NP[2]=2  Or  3:  As  an  alternative,  EQN  will  generate  a  complete  set  of  Hugoniot 
data  from  a  linear  or  quadratic  fit  for  either  Us-Up  (NPr21=21  or  P-Up  (NP[2)=31.  Either 
of  these  fits,  again,  may  be  comprised  of  up  to  10  segments.  The  initial  density, 

SETUP  [1],  and  the  three  coefficients  of  the  fit,  SETUP[2-4],  are  input  along  with  the 
ends  of  the  fit  and  the  geometric  increment  (AU)  in  Up.  (Defaults  are  Umin=.0001 , 

AU=1 .1 .  The  code  tests  on  PMX  so  if  Umax  is  not  input,  it  defaults  to  PMX=50  kb.) 

NPr21=4:  Because  the  shock  data  appear  to  be  consistent  with  a  piecewise 
'  linear  Us-Up  relationship,  many  experimenters  have  generated  such  fits  and  several 
xEOS  models  use  these  Us-Up  fits  to  obtain  the  Hugoniot.  CRT  has  built  a  data  base  of 
the  fits  that  can  be  used  to  generate  data  sets  for  comparison  with  other  models  or 
data.  Setting  NP[2]=4  will  cause  EQN  to  search  the  Us-Up  data  base  for  the  fit 
corresponding  to  an  ID  number  (ENUM=SETUP[1])_  Author  (IAUT=SETUP[2])  and 
initial  density  (RHOO=SETUP[3]).  After  finding  the  coefficients  of  the  Us-Up  fit 
requested,  EQN  generates  the  complete  Hugoniot  data  set  as  above. 


I 
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4.2.3  EXPHUG. 


For  ROUT=EXP,  subroutine  EXPHUG  is  selected  to  process  experimental 
Hugoniot  and  release  data  from  the  DNA  EOS  data  base.  Input  to  EXPHUG  is 
summarized  as  follows; 

NE 

1  LINE  Reference  number  of  data  for  plotting  or  other  routines 

2  LEX  if  *0;  generate  a  subset  of  data  from  a  previous  data  set  denoted 

by  NP(1) 

3  iHEL  if  *  0;  ignore  Hugoniot  Elastic  Limit  (HEL)  input 

1 ;  use  HEL  input  where  available  for  each  author 
2;  use  HEL  input  for  ALL  subsequent  data 

SETUE 

1  EDIT1  Number  of  experimental  Hugoniot  data  set 

if  EDIT1  >  100,  data  on  Unit  10  (Rocks.grout.H2O) 
if  "  <1 00,  data  on  Unit  1 1  (Metals,  elements) 

2  RHOMN  Minimum  initial  p  to  be  included  in  data  set  {0.} 

3  RHOMX  Maximum  initial  p  to  be  included  in  data  set  {1 00.} 

For  NPf21=0  HUGEXP  reads  data  from  the  files  containing  the  shock  and 
release  data  bases  and  calculates  the  complete  set  of  Hugoniot  parameters,  i.e.,  P,  p, 
Us,  Up,  and  E.  The  routine  searches  the  appropriate  data  base  corresponding  to  the 
'  material  number  EDIT1  {SETUP[1]).  Only  points  with  initial  densities  between  input 
*  limits  RHOMN  and  RHOMX  (SETUP[2  and  3])  are  stored  for  future  plotting.  If  the  limits 
are  not  specified,  all  the  data  for  the  selected  material  are  included. 

The  principal  Hugoniot  for  most  materials  is  not  a  single  shock  wave  over  the 
entire  range  of  pressures  of  interest.  The  stress  at  which  the  material  can  no  longer 
support  the  elastic  deviators,  solid-solid  or  solid-fluid  phase  changes  and  other 
dynamic  processes  (e.g.,  pore  collapse)  may  cause  a  decomposition  of  the  Hugoniot 
into  several  waves,  a  precusor  whose  pressure  is  refered  to  as  the  Hugoniot  elastic 
limit  (HEL)  and  the  main  shock3.  Since  the  precursor  may  be  due  to  a  phase  change  or 
other  process,  HEL  may  be  a  misnomer,  however,  it  is  used  to  describe  the  precursor 
stress  regardless  of  its  cause.  The  presursor  is  distingushed  by  the  fact  that  it  travels 


3  It  is  even  possible  to  decompose  a  shock  into  more  than  2  waves.  The  waveforms  in  the  marble  in  the  DISTANT 
MOUNTAIN  tests  suggest  a  3  wave  structure.  Such  a  wave  set  is  not  considered  here. 
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with  a  wave  speed,  Usp,  that  is  faster  than  that  of  the  main  wave.  If  a  precursor  is 
present,  the  main  wave  is  actually  impinging  on  preshocked  material  and  the  Hugoniot 
jump  conditions  must  be  altered  to  account  for  the  change.  Where  the  precursor  was 
reported,  it  is  included  in  the  data  base. 

One  problem,  particularly  in  older  data  obtained  in  less  sophisticated 
experiments,  is  that  the  precursor  may  have  been  overlooked.  In  consideration  of  this 
possiblity,  an  option  is  available  in  processing  the  data.  If  NP[3]=2  and  the  data  base 
includes  a  line  with  the  HEL  Usp,Upp  data,  these  data  are  used  to  determine  if  any 
subsequent  shock  points  are  multiple  shocks  (Us<Usp).  If  so,  the  complete  Hugoniot 
state  is  calculated  using  the  HEL  data  as  the  initial  state.  If  NP[3]=1 ,  only  the  data  from 
the  author  reporting  the  HEL  are  checked  and  used  where  appropriate.  If  NP(3)=0, 
only  experimental  data  for  which  the  precursor  was  reported  by  the  author  are  treated 
as  double  shocks. 

If  NPf2j>0.  EXPHUG  assumes  that  a  previous  call  has  stored  a  data  set  for  the 
material  in  line=NP[2]  and  this  call  is  to  generate  a  subset  of  that  data  limited  by 
RHOMX  and  RHOMX.  The  NP[2]>0  option  is  useful  in  that  it  preserves  the  symbol 
types  associated  with  each  experimenter.  This  provides  some  continuity  when  plotting 
the  various  data  sets. 

An  example  of  the  output  generated  by  ROUT=  EXP  is  presented  in  Table  4-3. 
The  initial  printout  lists  the  authors  and  the  Hugoniot  data  in  the  order  stored  in  the 
ROCKHUG  (or  METHUG)  database.  The  two  letter  designation  at  the  left  of  each  row 
designates  the  two  measured  variables,  e.g.,  UP  for  particle  velocity  and  pressure;  SP 
for  shock  velocity  and  pressure..  The  remainder  of  the  line  contains  the  complete  set 
of  data  associated  with  the  Hugonot  point.  If  a  release  path  for  the  data  point  is  also 
stored  in  the  RELAD  database,  a  second  line  denoting  that  is  printed  next  After  culling 
the  Hugoniot  data  in  the  database  and  eliminating  those  outside  the  requested  density 
limits,  the  code  sorts  the  remaining  points  and  produces  a  second  list  of  the  results  in 
order  of  increasing  pressure. 
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Table  4-3.  Example  of  the  output  generated  by  ROUT=EQN. 


Part  1  -  Echo  of  the  input  from  ROCKHUG 


PHYLLITE 
TYPE  ETA  P 

(MBAR) 


Us  Up 
(C/US)  (C/US) 


RHO 


init  RHO  ENG  1+PV/E  IDN- 
(TE/G) 


FURNISH  (93) 

UP  1.0800  0.0523  0. 

RA  RELEASE  PATH 
UP  1.1435  0.0861  0. 

RA  RELEASE  PATH 
UP  1.2170  0.1627  0. 

RA  RELEASE  PATH 
UP  1.4015  0.2658  0. 

RA  RELEASE  PATH 
UP  1.6847  0.6571  0. 

RA  RELEASE  PATH 
UP  1.8710  1.3971  1. 

RA  RELEASE  PATH 
KT  SLATE  (93) 

UP  1.0091  0.0044 


UP  1.0187 
UP  1.0332 
UP  1.0567 
UP  1.0798 
UP  1.0615 


0.0098 

0.0202 

0.0345 

0.0486 

0.0406 


5076 

FROM 

5004 

FROM 

5771 

FROM 

5819 

FROM 

7682 

FROM 

0466 

FROM 

4226 

4416 

4787 

4843 

4900 

,5057 


0.0376  2.959 
FILE  FOR  SHOT 
0.0628  3.133 
FILE  FOR  SHOT 
0.1029  3.335 
FILE  FOR  SHOT 
0.1667  3.840 
FILE  FOR  SHOT 
0.3122  4.616 
FILE  FOR  SHOT 
0.4872  5.126 
FILE  FOR  SHOT 


0.338 

SLP1 

0.319 

SLP2 

0.300 

SLP3 

0.260 

SLP4 

0.217 

SLP7 

0.195 

SLP6 


2.765  0.362 
2.791  0.358 
2.831  0.353 
2.895  0.345 
2.959  0.338 
2.909  0.344 


0.0038 
0.0081 
0.0154 
0.0260 
0.0362 
0.0293 

lines  missing 


2.740  7 
2.740  1 
2.740  5 
2.740  1 
2.740  4 
2.740  1 


2.740  7 
2.740  3 
2.740  1 
2.740  3 
2.740  6 
2.740  4 


069E-04 

1 

.  972E-03 
2 

.  294E-03 

3 

.  389E-02 

4 

.  873E-02 

5 

.  187E-01 

6 

.220E-06 
.  281E-05 
.  186E-04 
.380E-04 
.  552E-04 
.  292E-04 


FURNISH (92d)  froz 
UP  1.0215  0.0133  0.4799 

UP  1.0347  0.0254  0.5254 

UP  1.0540  0.0354  0.5012 

UP  1.0829  0.0495  0.4845 

UP  1.1150  0.0819  0.5379 

UP  1.0378  0.0239  0.4888 

UP  1.0736  0.0483  0.5073 

KTECH ( 92 )  t>0 
SP  1.0119  0.0084  0.5060 

SP  1.0222  0.0164  0.5200 

SP  1.0370  0.0268  0.5200 

SP  1.0519  0.0364  0.5140 

SP  1.0729  0.0530  0.5270 


0.0101 

0.0176 

0.0257 

0.0371 

0.0555 

0.0178 

0.0348 

0.0059 

0.0113 

0.0185 

0.0254 

0.0358 


2.797  0.358 

2.844  0.352 
2.894  0.345 
2.980  0.336 
3.057  0.327 

2.845  0.352 
2.939  0.340 

2.833  0.353 
2.852  0.351 
2.883  0.347 
2.935  0.341 
3.015  0.332 


2.738  5. 101E-05 
2.749  1.549E-04 
2.746  3.302E-04 
2.752  6. 882E-04 
2.742  1 . 540E-03 
2.741  1 . 584E-04 
2.737  6 . 055E-04 


26.003  13- 
1 

14.935  13- 
2 

10.216  13- 

3 

5.982  13- 

4 

3.921  13- 

5 

3.296  13- 

6 

0.000  14- 
0.000  14- 
61.171  14- 
36.252  14- 
26.071  14- 
33.520  14- 


94.022  16- 
58.704  16- 
38.003  16- 
25.120  16- 
18.384  16- 
53.925  16- 
28.156  16- 


2.800  1 . 758E-05  0.000  17- 
2.790  6 . 389E-05  91.002  17- 
2.780  1 . 718E-04  55.098  17- 
2.790  3 . 221E-04  39.500  17- 
2.810  6 . 405E-04  28.450  17- 


lines  missing 


Part  2  ~  Hugoniot  data  arranged  by  increasing  PH 


PHYLLIT  (2.72  -  2.82) 


N 

IDN 

ETA 

P 

Us  Up 

RHO 

V 

(MBAR) 

(C/US)  (C/US) 

1 

7 

14 

1.0091 

0.0044 

0.4226  0.0038 

2.765 

0.362 

2 

35 

18 

1.0118 

0.0083 

0.5050  0.0059 

2.823 

0.354 

3 

30 

17 

1.0119 

0.0084 

0.5060  0.0059 

2.833 

0.353 

4 

38 

18 

1.0119 

0.0085 

0.5070  0.0059 

2.853 

0.350 

5 

8 

14 

1.0187 

0.0098 

0.4416  0.0081 

2.791 

0.358 

6 

16 

15 

1.0146 

0.0109 

0.5283  0.0076 

2.762 

0.362 

lines  missing 

37 

20 

15 

1.1168 

0.0865 

0.5490  0.0574 

3.066 

0.326 

38 

15 

14 

1.1425 

0.1028 

0.5485  0.0684 

3.130 

0.319 

39 

3 

13 

1.2170 

0.1627 

0.5771  0.1029 

3.335 

0.300 

40 

4 

13 

1.4015 

0.2658 

0.5819  0.1667 

3.840 

0.260 

41 

5 

13 

1.6847 

0.6571 

0.7682  0.3122 

4.616 

0.217 

42 

6 

13 

1.8710 

1.3971 

1.0466  0.4872 

5.126 

0.195 

for 

ROUT 

=  EXP  User  LINE 

*  2  equal  to 

driver 

Lin  2 

RHOI  ENG 

(TE/G) 

2.740  7 . 220E-06 
2.790  1 . 735E-05 
2.800  1 . 758E-05 
2.820  1.767E-05 
2.740  3. 281E-05 
2.722  2 . 888E-05 

2.745  1 . 647E-03 
2.740  2 . 339E-03 
2.740  5 . 294E-03 
2.740  1 . 389E-02 
2.740  4.873E-02 
2.740  1 . 187E-01 


N 


1  SLP-1 

2  SLP-2 

3  SLP-3 

4  SLP-4 

5  SLP-7 

6  SLP-6 


7  3621 

8  3616 

9  3617 

10  3618 

11  3619 

12  3622 


23  PFP1A 

24  PFP2 

25  PFP3 

26  PFP4 

27  PFP5 

28  PFR1 

29  PFR2 

30  3602 

31  3586 

32  3584 

33  3585 

34  3587 


1+PV/E 


0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

18.128 

15.038 

10.216 

5.982 

3.921 

3.296 
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4.2.4  FLIP  -  Reverse  Hugoniots. 

When  a  shock  wave  passes  through  two  different  materials  or  when  one  material 
strikes  the  other,  it  is  possible  to  determine  the  shock  conditions  at  the  interface 
knowing  only  one  parameter  if  the  Hugoniots  of  the  two  media  are  also  known.  This  is 
done  by  rotating  one  of  the  Hugoniots,  Ph  about  the  Up=0  axis  and  than  translating  the 
rotated  (flipped)  curve,  PF  so  that  the  PF=0  point  remains  on  the  x-axis.  Calling  FLIP 
provides  the  tools  to  flip  and  translate  the  Hugoniots  for  several  types  of  input,  namely. 


LINE  Reference  number  of  data  for  plotting  or  other  routines 
LI  Previous  line  to  be  flipped 
L2  Previous  line  to  compare  to  LI  to  obtain  BCF's 
NPP  read  in  additional  Up,P  data  for  flipped  Hugoniots 


For  NP[3]  =  0 

0;  SETUP  =  up  to  12  Up's  where  Ph  and  PF  intersect 

NP[4]  =1  .  P's  "  "  "  "  " 

2  "  "  6  pairs  of  Up,P  that  PF's  pass  through 


SETUP  NP[3]>0 

1  P2MAX  Max  pressure  on  L2  for  BCF's  [def=. 25  Mb] 

2  PINC  Ratio  of  successive  pressures  along  L2  for  BCF's 


For  NPr31=0  FLIP  flips  the  data  set  defined  by  NP[2]  and  translates  the  result  to 
intersect  with  the  original  curve.  If  NP[4]  =  0  or  1 ,  the  Hugoniot  is  translated  so  that  PF 
'intersects  PH  at  either  Up  (NP[4]=0)  or  P  (NP[4]=1),  where  Up=SETUP[i],i=1,12  or 
P=SETUP[i],i=1,12.  If  NP[4]=2,  a  set  of  flipped  Hugoniots  are  generated  that  intersect 
with  up  to  six  Up-  P  pairs  read  into  SETUP[1-12],  Additional  pairs  are  again  read  by 
setting  NP[8]*0.  For  each  line  flipped,  a  line  of  text  is  printed  given  the  point  of 
rotation,  thusly; 


ME191  r  =  1.91 

FOR  LINE  11  HUGONIOT  IS  FLIPPED  ABOUT  Up=  0.100  (cm/us) 
User  line  6  copy  of  line  4  and  flipped  to  make  6  new  lines 


,  For  NPf31>0.  FLIP  will  generate  the  borei lole  correction  factors  relating  one 
medium  to  another.  The  BCF  is  defined  as  the  ratio  of  the  Up  or  P  in  the  second 
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material  to  the  corresponding  value  in  the  first.  The  shock  is  assumed  to  go  from  the 
data  set  defined  by  LI  {NP[2]}  into  the  one  defined  by  L2  {NP[3]>.  If  the  impedance 
(l=p0Us2)  of  both  media  are  identical,  the  BCFs  are  1.0.  If  the  impedance  is  less  in  the 
second  medium  [l2<li],  the  BCFjj  is  greater  than  1 .0  and  the  BCFp  is  less  than  1 .0.  If 
the  second  medium  has  a  higher  impedance  [l2>li],  the  BCFy  is  less  than  1.0  and  the 
BCFp  is  greater  than  1 .0. 


FLIP  finds  the  BCFs  starting  at  PTST=1  kbar  and  increasing  the  pressure  by 
10%  each  step.  PTST  and  the  corresponding  U  on  the  second  Hugoniot  are  one  pair 
of  data.  The  Hugoniot  of  the  first  medium  is  translated  until  it  intersects  the  second  at 
PTST  and  the  pressure  and  particle  velocity  at  its  intersection  with  the  original  first 
medium's  Hugoniot  are  the  other  two  values  needed  to  determine  the  ratios.  A  printout 
for  each  PTST  is  generated  with  the  two  pressures  and  BCFp  and  the  two  velocities 
with  BCFu  up  to  a  maximum  PTST  of  P2MAX  [def=250  kbars],  Table  4-4. 

Table  4-4.  Example  of  FLIP  printout  for  borehole  correction  factors. 

BOREHOLE  CORRECTION  FACTORS  between  ME191  and  HTH3 


N1 

N2 

Pff 

Pbh 

bcfP 

Uff 

Ubh 

bcfUp 

44 

31 

0.0009 

0.0010 

0.8415 

0.0033 

0.0029 

1.1569 

45 

33 

0.0010 

0.0012 

0.8503 

0.0037 

0.0033 

1.1483 

46 

34 

0.0011 

0.0013 

0.8535 

0.0040 

0.0035 

1.1441 

47 

35 

0.0012 

0.0014 

0.8567 

0.0042 

0.0037 

1.1401 

48 

36 

0.0013 

0.0015 

0.8598 

0.0045 

0.0040 

1.1362 

50 

38 

0.0014 

0.0017 

0.8625 

0.0049 

0.0044 

1.1308 

51 

39 

0  .‘0015 

0.0018 

0.8638 

0.0052 

0.0046 

1.1281 

84 

70 

'0.0180 

0.0204 

0.8827 

0.0354 

0.0321 

1.1060 

92 

79 

0.0362 

0.0401 

0.9038 

0.0593 

0.0549 

1.0816 

122 

110 

0.1538 

0.1599 

0.9618 

0.1639 

0.1595 

1.0274 

124 

112 

0.1679 

0.1737 

0.9668 

0.1740 

0.1700 

1.0235 

127 

115 

0.1885 

0.1934 

0.9743 

0.1884 

0.1848 

1.0196 

130 

117 

0.2049 

0.2100 

0.9755 

0.2003 

0.1967 

1.0183 

132 

119 

0.2244 

0.2296 

0.9773 

0.2139 

0.2103 

1.0168 

I 
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4.2.5  HUG  -  The  Hugoniot  Driver. 

For  ROUT=HUG  ,  subroutine  HUG  is  selected  to  generate  the  principal 
Hugoniot4  for  the  material  specified.  Input  to  HUG  is  limited  to  several  optional 
parameters,  namely  for  ROUT  =  HUG 

NE 

1  LINE  Reference  number  of  data  for  plotting  or  other  routines 

SETUP 

1  PUP  Max  stress  to  load  Hugoniot;  {100Mb}5 

2  YVEL  Initial  velocity  (volume)  step  size;  {dv  =  0.001  cm/jis} 

3  GSY  Initial  Hydrostatic  pressure  {0  Mb} 

HUG  begins  by  calling  EQSINIT  to  generate  a  unit  cell  compatible  with  the 
CRALE2  code  subroutines.  Using  the  CRALE  routines  insures  consistency  between 
the  results  of  the  driver  and  the  2D  code  calculations.  The  Hugoniot  is  found  by  setting 
the  vertical  velocities  at  the  top  of  this  cell  to  a  negative  value  thus  generating  a 
uniaxial  compressive  strain  condition.  (Defaults  of  At=1  and  YVEL=. 001,  result  in  an 
initial  0.1%  change  in  volume.  The  time  step  is  just  used  to  generate  displacements 
and  does  not  imply  a  strain  rate!)  The  data  generated  by  HUG  are  tagged  with  a 
number  between  1  and  99,  NP(1),  so  that  they  can  be  retrieved  for  use  in  plots  or  other 
,  routines,  e.g.,  EQN  and  FLIP. 

\ 

To  obtain  the  complete  Hugoniot,  HUG  runs  through  a  series  of  steps  for 
continually  increasing  compressions.  Each  step  begins  by  adjusting  YVEL  so  that  the 
next  Hugoniot  stress,  oy,  is  at  least  1.05  but  less  than  1.25  times  the  previous.  YVEL 
may  also  be  adjusted  if  the  code  had  trouble  converging  on  the  last  answer.  HUG 
converges  on  the  Hugoniot  stress  by  fixing  the  change  in  strain,  ey,  and  iterating  on  the 
change  in  the  cell's  specific  energy  until  it  differs  by  less  than  0.01%  from  the  energy 
density  in  a  strong  shock  as  determined  by  the  Hugoniot  jump  condition,  i.e. , 

EH=.5aydey  (4.2) 


4  The  principal  Hugoniot  is  the  locus  of  shock  states  for  a  material  whose  pre-shock  state  is  P=0  Un=0  E=E„  and 
'  P=P0 

0  {}  indicates  default  value 
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where  EH  is  the  energy/gram;  cy,  the  axial  (vertical)  stress;  and  ey,  the  axial  strain. 

Since  the  motions  are  uniaxial;  ex  =  e*  =  exy  =  0,  so  the  compression  and  axial  strain 

are  related  by 

p  =  -ey/(1.-ey)  (4-3) 

and  Equation  4.1  is  equivalent  to  the  more  familiar  but  less  general 

Eh=.5PAV  (4-4) 

for  a  hydrodynamic  material.  Having  converged  on  the  Hugoniot  pressure  and  energy 
for  the  current  strain,  ey,  HUG  stores  and  prints  the  shock  state  parameters  of  interest 
(p,  a,  PH,  K,  G,  Us,  Up,  etc.)  plus  the  pressure  corresponding  to  p  and  E=E0,  i.e.,  the 
zero-energy  hydrostat.  The  routine  proceeds  to  the  next  strain  until  the  Hugoniot 
stress  exceeds  PUP  [default  =  100  Mb,  can  be  set  by  SETUP(2)]. 

Upon  convergence  at  each  step,  HUG  checks  that  the  shock  velocity,  Us,  is  still 
a  maximum.  If  Us  is  less  than  that  attained  at  a  lower  stress,  the  shock  will  not 
propagate  as  a  single  wave,  but  breaks  into  a  precursor,  traveling  at  the  previous 
higher  velocity  and  a  second  front  moving  with  a  lower  speed.  The  Hugoniot  equations 
for  such  a  double  shock  are  slightly  more  complex  since  the  initial  conditions  for  the 
second  pulse  is  the  state  (i.e.,  P,  p,  Up,  E)  of  the  material  behind  the  precursor.  When 
HUG  finds  s.uch  a  condition,  the  routine  saves  the  data  at  the  first  front,  i.e.,  the 

r 

Hugoniot  Elastic  Limit  (HEL)  and  Equation  4.2  is  modified  appropriately.  Note,  the  drop 
'in  wave  speed  could  result  from  a  phase  change,  void  collapse  or  other  non-elastic 
process  as  well  as  the  onset  of  plastic  yielding,  so  the  term  HEL  may  be  a  misnomer  for 
some  materials.  Nonetheless  HEL  is  used  to  label  the  initial  point  of  the  double  shock 
region  regardless  of  the  cause.  HUG  continues  to  compare  the  current  shock  velocity, 
Us,  against  the  maximum  and  when  Us  equals  or  exceeds  the  previous  high,  the  shock 
is  again  assumed  to  merge  into  a  single  wave  starting  from  the  original  state  of  the 
material.  While  even  more  complex  shock  states  are  possible,  HUG  only  considers  at 
most  a  two-wave  structure. 

If  the  EOS  being  exercised  is  very  sensitive  to  changes  in  energy,  it  is  possible 
for  the  Hugoniot  to  reach  a  maximum  compression  or  even  double-back,  in  P-V  or  P-|i 
space.  Gases  and  highly  hysteretic  solids  such  as  snow  texhibit  this  behavior.  In  this 
case  the  test  compression  may  not  be  attainable  in  a  single  shock  and  the  standard 
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(energy)  iteration  can  not  converge.  When  HUG  cannot  converge  within  30  iterations 
or  if  the  energy  changes  from  the  last  step  by  more  than  a  factor  of  2,  the  code 
proceeds  by  guessing  a  change  in  E  and  iterates  on  density  for  an  additional  20  tries. 

If  it  still  cannot  converge,  the  last  five  steps  are  printed  out  and  the  code  continues. 
HUG  assumes  the  material  is  passing  through  some  unique  local  phase.  While  the 
routine  may  not  be  able  to  converge  in  a  small  range  of  densities  due  to  local 
phenomena,  usually  it  will  converge  again  at  higher  densities  and  energies. 

After  each  step,  HUG  automatically  generates  a  one  line  printout.  Output  from  a 
typical  run  is  shown  at  the  end  of  this  subsection.  The  printout  for  each  step  includes 
as  much  information  about  the  state  of  the  material  as  can  be  squeezed  onto  a  single 
132  character  line,  beginning  with  p,  p,  the  pressure  at  E=Eo,  Ph,  oy  and  the  strain,  ey. 
The  printout  also  includes  the  various  wave  speeds  [shock  (Us),  local  compressional  or 
sound  (Cp),  and  shear  (Cs)],  the  bulk  (K)  and  shear  (G)  moduli,  Poisson's  ratio  (v),  the 
ratio  of  Vo/ Y  (DSIG),  the  yield  surface  (Y),  the  plastic  strain  (EPL,  epi)  and  Hugoniot 
energy,  Eh.  When  a  material  melts  (i.e.,  E  >  ES),  G,  Y  and  ep|  are  set  to  zero  and  an 
effective  gamma  (GAM),  defined  as 

Y  =  l.+P/pE  (4.5) 

and  the  specific  energy  density  (EH)  are  printed  instead  of  v  and  ep|. 


t 


l 


46 


w 

o 

X 


a 

x 

\ 

O' 

C 

La 


o 

o 


X 

a: 

o 

DC 

a: 


o 

o 


O 

X 

X 


CO 

Q 


D 

2 


I 
I 
I 

0  0 


m 

2 


O  I 
i 
l 
I 
l 

a  i 

0 

CO 

ra 


a 

D 


O  CM 
O  -H 
+  I 

La  x 

o  cn 
O  tp 

o  m 

O  *H 
O  'C* 

o  o 

+  I 
X  X 

O  CM 

o  o 


o  o 
o  o 


o  o 

CM  CM 

m  m 

o  o 

in  in 

-<r  ^ 

o  o 

CM  CM 

o  o 

o  o 
o  o 
o  o 
m  in 

o  o 

m  m 

p  «h 

cm  in 


co  cn 
m  in 


o  cn 

o  in 

o  o 

o  o 
o  o 
o  o 
o  o 

d  d 

o  0 
o  o 
h-  i 
La  La 

O  tp 
-  O  0 
O  CM 

o  cn 

I 

O  ID 
CM  O 
I  t 
l:  La 
o  cn 
o  p 
o  in 


m 

in 

00 

tp 

in 

m 

m 

in 

00 

CD 

CD 

00 

00 

cn 

cn 

cn 

cn 

cn 

VD 

cn 

cn 

cn 

o 

o 

o 

o 

rH 

*H 

rH 

rH 

in 

m 

m 

in 

in 

in 

m 

m 

m 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

rH 

rH 

•H 

rH 

rH 

rH 

rH 

o 

o 

w—4 

rH 

o 

o 

o 

o 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

o 

o 

o 

o 

o 

i 

i 

i 

1 

i 

i 

i 

i 

i 

i 

i 

I 

i 

i 

i 

1 

1 

» 

i 

i 

1 

+ 

H- 

+ 

+ 

La 

x 

La 

La 

La 

La 

La 

La 

LO 

ua 

La 

La 

La 

ua 

La 

La 

La 

La 

La 

LO 

LO 

La 

LO 

La 

La 

La 

La 

m 

p 

VD 

CM 

o 

r~- 

m 

cn 

CO 

o 

CO 

TP 

rH 

cn 

VD 

rH 

■'T 

T 

VD 

cn 

o 

CD 

rH 

rH 

o 

HP 

00 

cn 

m 

CM 

TT 

CM 

CM 

in 

cn 

VD 

VD 

rH 

o 

cn 

r- 

r- 

CM 

cn 

CM 

CM 

rH 

cn 

m 

CD 

rH 

o 

Tp 

p 

o 

p 

CM 

CD 

o 

CM 

m 

00 

rH 

m 

cn 

cn 

cn 

VD 

O 

cn 

in 

CM 

VD 

cn 

CM 

VD 

CM 

CM 

CM 

cn 

cn 

TV 

m 

VD 

«c 

"T 

in 

m 

in 

VD 

m 

cn 

in 

VD 

p- 

cn 

rH 

rH 

cn 

cn 

m 

m 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

m 

in 

m 

m 

in 

in 

m 

in 

m 

m 

m 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

i 

LO 

t 

aa 

i 

x 

La 

La 

u 

U3 

u 

ill 

Ll) 

La 

LO 

La 

LO 

La 

La 

LO 

LO 

La 

LO 

La 

La 

La 

La 

La 

LO 

L0 

O 

Tp 

r* 

o 

CM 

TP 

in 

VD 

VO 

CM 

o 

r- 

Tp 

•H 

00 

r- 

r- 

r- 

r~ 

p* 

r- 

p- 

r- 

r~ 

r~ 

P 

m 

cn 

cn 

tt 

TP 

rH 

rH 

o 

o 

o 

cn 

cn 

r-~ 

r- 

r~ 

r~ 

p~ 

r~ 

r~ 

r* 

p~ 

p- 

P 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

m 

cn 

CM 

CM 

VD 

VD 

VD 

VD 

VD 

VD 

VO 

VO 

VD 

VD 

VD 

«H 

CD 

m 

m 

m 

in 

m 

m 

m 

m 

in 

in 

in 

m 

in 

m 

O 

o 

O 

o 

O 

o 

o 

o 

O 

O 

O 

o 

o 

rH 

H 

H 

rH 

*H 

H 

•H 

rH 

rH 

r- 1 

rH 

rH 

rH 

rH 

O 

o 

O 

o 

O 

o 

o 

o 

o 

O 

O 

rH 

•H 

H 

rH 

rH 

rH 

H 

rH 

rH 

rH 

*H 

rH 

rH 

rH 

rH 

rH 

O 

o 

O 

o 

o 

o 

o 

o 

o 

O 

O 

rH 

CM 

CM 

CM 

CM 

CM 

o 

rH 

cn 

in 

VD 

GO 

rH 

O 

o 

O 

o 

o 

o 

o 

o 

o 

o 

o 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

m 

m 

m 

in 

in 

m 

VD 

O 

o 

O 

o 

o 

o 

o 

o 

o 

o 

o 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

m 

m 

m 

in 

m 

m 

in 

m 

m 

in 

m 

in 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

p 

cn 

o 

00 

m 

m 

rH 

p 

tp 

cn 

cn 

r- 

CD 

m 

cn 

o 

o 

r- 

co 

CD 

cn 

m 

.-H 

HP 

rH 

cn 

m 

m 

VD 

p 

CD 

00 

cn 

cn 

cn 

cn 

CM 

rH 

o 

co 

cn 

rH 

rH 

cn 

in 

o 

HP 

CO 

cn 

VD 

o 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

m 

m 

m 

in 

m 

m 

TT 

in 

rH 

VD 

o 

in 

o 

VD 

CM 

cn 

p 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

r* 

r* 

VD 

VD 

in 

in 

r- 

p 

VD 

VD 

O 

O 

o 

O 

o 

O 

O 

O 

o 

o 

o 

o 

o 

o 

o 

o 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

rH 

rH 

rH 

rH 

o 

o 

o 

CM 

tp 

00 

cn 

O 

00 

m 

o 

cn 

cn 

cn 

cn 

r* 

TT 

rH 

r~ 

in 

O 

cn 

cn 

rH 

r> 

P 

VD 

VD 

00 

o 

CM 

TP 

VD 

0 

rH 

CM 

o 

CM 

cn 

m 

CO 

o 

<n 

CD 

CM 

m 

VD 

rH 

CD 

HT 

CM 

m 

r- 

CM 

CM 

CM 

cn 

cn 

cn 

cn 

cn 

Tp 

TP 

VD 

p 

00 

cn 

o 

CM 

cn 

m 

r* 

HP 

cn 

VD 

co 

m 

VD 

m 

in 

m 

m 

m 

m 

m 

in 

m 

o 

o 

o 

o 

rH 

rH 

rH 

rH 

cn 

in 

r* 

o 

cn 

VD 

CM 

cn 

cn 

o 

o 

o 

o 

o 

o 

o 

o 

o 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

*T 

m 

in 

in 

CD 

HP 

rH 

cn 

CM 

cn 

HP 

p 

o> 

rH 

00 

rH 

CO 

o 

p 

00 

p 

rH 

00 

VD 

CD 

r- 

o 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

CM 

CM 

cn 

CM 

TP 

cn 

VD 

in 

CO 

o 

CM 

cn 

in 

r* 

as 

o 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

m 

m 

m 

in 

in 

m 

in 

*H 

«H 

o 

rH 

*H 

rH 

rH 

rH 

•H 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

O 

O 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

p 

CD 

0 

o 

rH 

CM 

cn 

tp 

m 

CM 

m 

p 

O 

cn 

VD 

00 

cn 

VD 

m 

CD 

•a- 

m 

o 

rH 

p»* 

p 

CD 

tp 

tp 

HT 

m 

in 

in 

in 

m 

m 

in 

m 

m 

VD 

VD 

VD 

VD 

<n 

VD 

cn 

CM 

VD 

o 

in 

HP 

in 

o 

VD 

m 

m 

m 

m 

m 

in 

in 

m 

m 

p 

p 

r~~ 

p 

p 

r- 

r- 

CM 

CM 

CM 

cn 

cn 

CM 

m 

cn 

CM 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

cn 

n 

cn 

cn 

«H 

cn 

VD 

p 

CM 

TP 

CO 

o 

00 

CM 

CO 

p 

cn 

cn 

cn 

'O’ 

<0* 

rH 

o 

CM 

CM 

CO 

V 

VO 

cn 

P 

m 

m 

VD 

m 

tp 

cn 

VD 

cn 

cn 

VD 

CM 

cn 

VD 

cn 

rH 

cn 

O 

CM 

o 

in 

cn 

o 

r- 

o 

rH 

VD 

tp 

Tp 

tt 

VO 

VD 

VD 

VO 

CM 

cn 

rH 

CM 

cn 

cn 

TP 

m 

VD 

-<r 

cn 

o 

VD 

cn 

CM 

VO 

CD 

in 

in 

in 

TP 

tp 

TP 

tp 

tp 

cn 

in 

in 

m 

lO 

in 

m 

in 

m 

m 

VD 

r* 

P* 

00 

cn 

rH 

CD 

VD 

o 

o 

o 

o 

O 

o 

O 

O 

o 

o 

o 

o 

o 

o 

o 

o 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

m 

VD 

VD 

P 

TP 

0 

m 

0 

TP 

0 

m 

cn 

CM 

o 

cn 

cn 

CM 

VD 

CM 

*H 

00 

cn 

P* 

cn 

CD 

r- 

rH 

VD 

O 

CM 

P 

VD 

VD 

p 

p 

00 

CO 

0 

cn 

o 

CM 

nr 

p 

rH 

TP 

CO 

CM 

<n 

VD 

VD 

rH 

CM 

o 

in 

cn 

■*r 

rH 

CM 

O 

o 

o 

o 

o 

o 

o 

O 

rH 

cn 

cn 

cn 

O 

o 

o 

rH 

rH 

in 

cn 

cn 

m 

rH 

00 

HP 

in 

P 

O 

o 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

rH 

rH 

rH 

rH 

00 

00 

CD 

cn 

cn 

o 

rH 

cn 

m 

rH 

P 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

d 

o  O 

rH 

rH 

cn 

m 

in 

CM  CM  CM  CM 
O  O  O  O 
till 

La  La  La  La 
cn  co  cn  m 
cn  r-  0  0 

o  h  cm  m 


CM 

o 

I 

La 

CM  ID 
O  ID 


CM  CM  CM  CM  CM 

O  O  O  O  O 

1  1  I  I  I 

La  La  La  La  La 

CM  «D  CD  H  m 

r-  as  o  cn  id 

in  id  ao  o')  o 


La 
ID 

VD  I  I  I  I  I 


o  o  o  o 

I  I  l  I 
La  La  La  La 
0  cn  0  m 
vd  cn  m  p 
r-  r*  cd  cd 


o  o  o 

•Hill 

o  La  La  La 

I  CM  ?  CD 

La  h  in  ct) 
0  0  0  0 
CD  •  •  • 

H  H  H 


»H  *H  H  *H  La 

o  o  o  o  o 

I  I  I  I  VD 

La  La  La  La  p 

CM  V  H  H  • 

rH  0  0  TP  CM 

m  cn  in 


o  o  o  o  o  o  o 

i  i  i  i  i  i  i 

La  La  La  La  ua  La  La 

in  m  vo  vd  o  io  ci 

CM  H  O  00  O  CD  TT 

vd  p  CD  cm  tp  nr  lO 


in  in  in  in  o 


000 

i  i  i 


r-  r-  p 


p  tp  o  T~t 


h  cm  cm  n  n  'T  cncn^TTr 


> 

o 

£>  i 


•r  id  cn  n  c  cm  oo  c  in  vd  h 


iH 

o 

VD 

cn 

cn 

ro 

<n 

0 

CM 

o 

o 

o 

O 

o 

tr* 

l*H 

cc 

1 

1 

i 

i 

i 

i 

t-H 

0 

< 

La 

LO 

La 

LO 

LO 

ua 

X 

0 

X 

m 

o 

m 

o 

00 

T 

co 

♦-H 

0 

X 

2 

o 

CD 

p 

CM 

CM 

HP 

X 

LO 

— 

o 

VD 

VD 

rH 

VD 

rH 

OC 

Oi 

• 

0 

E-h 

i 

cn 

rH 

in 

VD 

VD 

P 

t— t 

0 

1 

H 

i 

o 

VD 

cn 

cn 

m 

cn 

0 

i 

o 

O 

o 

o 

o 

o 

< 

< 

o 

i 

+ 

1 

i 

t 

i 

t 

La 

II 

i 

LO 

L0 

La 

LO 

LO 

La 

La 

X 

LO 

i 

O 

CD 

HP 

O 

cn 

o 

i 

o 

VO 

VD 

rH 

in 

rH 

E-* 

\ 

Cu 

i 

O 

i 

o 

rH 

m 

VD 

VD 

P 

w 

2 

O 

o 

m 

co 

P 

m 

rH 

O 

M 

o 

o 

HP 

• — l 

m 

HP 

0 

2 

o 

o 

o 

cn 

00 

CO 

X 

o 

o 

o 

rH 

rH 

CM 

cn 

X 

0 

o 

o 

rH 

rH 

D 

D 

o 

o 

O 

o 

o 

o 

(D 

X 

E 

• 

* 

• 

-C 

o 

o 

o 

o 

o 

o 

H 

o 

o 

o 

CM 

m 

CD 

rH 

rH 

HT 

HP 

HP 

HP 

O 

p 

P 

P 

P 

P 

P 

X 

X 

CM 

CM 

CM 

CM 

CM 

(M 

c 

rH 

CM 

p 

CD 

cn 

O 

p 

P 

p 

CD 

cn  cn  m  cn  cm 

o  o  o  o  o 

i  i  i  i  i 

La  La  La  La  La 

r-  o  o  co  rH 


o  o  o  o 

i  i  t  i 

La  La  La  La 


00 

CO 

cn 

cn 

rH 

rH 

rH 

rH 

rH 

cn 

cn 

cn 

cn 

cn 

rH 

rH 

rH 

rH 

o 

o 

o 

O 

o 

o 

o 

O 

o 

a  tla 

LO 

La 

LO 

LO 

X 

X 

X 

X 

cn 

m 

rH 

CM 

rH 

cn 

"<T 

cn 

CM 

CD 

in 

P 

cn 

rH 

rH 

CM 

CM 

CO 

CO 

cn 

<n 

cn 

rH 

rH 

rH 

o 

CD 

p 

CM 

VD 

in 

co 

cn 

o 

p 

HP 

co 

cn 

VD 

cn 

cn 

Ln 

tp 

m 

rH 

cn 

vd 

O 

cn 

cn 

p 

TP 

m 

p 

co 

cn 

rH 

cn 

HP 

o 

4  rH 

t-H 

rH 

rH 

CM 

rH 

rH 

CM 

cn 

o  o 

o 

o 

o 

O 

CM 

CM 

CM 

CM 

o  o 

o 

o 

o 

o 

o 

o 

O 

o 

o  »n 

VD 

o 

cn 

p 

rH 

in 

cn 

TP 

■>  m 

m 

VD 

VD 

VD 

cn 

o 

rH 

m 

-  p 

p 

p 

p 

P 

CM 

cn 

cn 

cn 

si  CM 

CM 

CM 

CM 

CM 

cn 

cn 

cn 

m 

H  CM 

co 

'p 

m 

m 

cn 

o 

rH 

CM 

0  00 

CD 

00 

co 

CD 

CM 

cn 

cn 

cn 

t-H 

rH 

rH 

rH 

cn  cn  cn  ■'y  (ovinininr-ciH 


o  o  o 

lit 

La  La  l: 


o  o 
o  o 


La  x 

in  0 


o  o  o  o  o  o 

o  o  o  o  o  o 

+  +  +  +  +  + 

x  x  x  x  x  x 

m  cd  cm  cd  r  in 

n  n  oo  cd  h  cm 


I: 

u  • 
o 

X  ' 


rH 

rH 

rH 

rH 

rH 

w 

(T3 

rH 

rH 

CM 

CM 

cn 

cn 

X 

0 

in 

cn 

CM 

TT 

in 

CL  rH 

o 

TP 

VD 

0 

0 

cn 

TP 

CM 

VD 

CM 

rH 

vD 

in 

rH 

Ln 

P 

TP 

p 

TP 

0 

VD 

cn 

m 

cn 

TT 

T5 

in 

VD 

rH 

0 

cn 

TP 

cn 

>  CM 

cn 

CM 

VD 

CM 

CM 

•H 

0 

CM 

TP 

TP 

in 

0 

o 

TP 

TP 

cn 

VD 

o 

TP 

f—H 

0 

cn 

0 

0 

TT 

P 

0 

si  CM 

CM 

rH 

rH 

CM 

CM 

0 

CM 

cn 

cn 

0 

0 

cn 

o 

a  o 

O 

rH 

rH 

rH 

rH 

w 

f-H 

i— i 

t-H 

CM 

CM 

CM 

<n 

H  0 

VD 

o 

rH 

o 

P 

0 

TP 

rH 

rH 

in 

cn 

0 

0 

a  v.o 

0 

0 

P 

p 

P 

cn 

CM 

VD 

0 

CM 

P 

m 

a  cn 

cn 

p 

0 

cn 

O 

c 

rH 

m 

TP 

cn 

TP 

P 

o 

■>  cn 

cn 

m 

m 

m 

VD 

C7 

VD 

VD 

cn 

O 

o 

t-H 

3 

rH 

rH 

a  tp 

m 

cn 

cn 

o 

O 

X 

rH 

rH 

CM 

TP 

m 

0 

P 

i  cn 

cn 

0 

0 

cn 

cn 

cn 

cn 

cn 

o 

o 

o 

o 

H  rH 

t-H 

t-H 

rH 

rH 

rH 

rH 

rH 

rH 

CM 

CM 

CM 

CM 

c 

•H 


cn 

cn 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

rH 

rH 

rH 

rH 

0 

rH 

rH 

o 

o 

o 

O 

X 

o 

o 

o 

rH 

rH 

rH 

CM 

M 

o 

o 

o 

o 

o 

o 

O 

O 

o 

O 

O 

o 

o 

O 

o 

o 

o 

o 

o 

o 

o 

o 

cn 

o 

o 

o 

O 

O 

O 

O 

G> 

i 

i 

1 

1 

rH 

1 

1 

1 

i 

| 

1 

1 

1 

1 

1 

rH 

1 

1 

i 

4- 

+ 

+ 

H- 

0 

4- 

H- 

+ 

+ 

+ 

H- 

+ 

> 

X 

X 

X 

X 

• 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

O 

X 

X 

X 

X 

X 

X 

X 

-H 

o 

0 

CM 

cn 

rH 

CM 

TP 

cn 

0 

o 

rH 

rH 

0 

0 

rH 

o 

0 

0 

Cn 

rH 

o 

cn 

cn 

p 

0 

CM 

CM 

P 

0 

u 

m 

m 

CM 

o 

0 

CM 

0 

m 

0 

O 

o 

m 

rH 

0 

a? 

0 

CM 

o 

CM 

cn 

cn 

cn 

cn 

TP 

cn 

o 

O 

0 

o 

cn 

V 

H  rH 

rH 

cn 

cn 

cn 

TP 

0 

X 

TP 

0 

0 

0 

p 

cn 

rH 

to 

H  *-H 

rH 

o 

o 

o 

o 

o 

o 

o 

rH 

rH 

»H 

CM 

rH 

>  o 

o 

o 

o 

o 

o 

* 

o 

o 

o 

o 

o 

o 

O 

co 

i 

i 

+ 

+ 

+ 

H- 

X 

+ 

H- 

H- 

+ 

+ 

+ 

H- 

3 

a  x 

X 

X 

X 

X 

X 

0 

X 

X 

X 

X 

X 

X 

X 

cr 

>  o 

rH 

cn 

rH 

o 

cn 

r~ 

cn 

p 

0 

CM 

CM 

p 

0 

a> 

4  cn 

P 

CM 

cn 

cn 

cn 

U 

Tp 

cn 

o 

o 

0 

o 

0 

)  0 

0 

TP 

0 

cn 

cn 

P 

CM 

0 

cn 

0 

0 

rH 

m 

La 

2 


0 

3> 

X 


H 

D 

O 

oc 


47 


CPU  TIME  USED  2.46815 


4.2.6  REL  --  Release  Adiabats  from  the  Hugoniot. 

After  a  call  to  HUG  has  produced  a  Hugoniot,  setting  ROUT  =  REL  will  call 
subroutine  REL  to  generate  the  paths  (release  adiabats)  followed  during  unloading. 
The  input  to  REL  are  summarized  as  follows; 

NE 

1  -  not  used  by  REL 

2  =  0  -  generate  release  adiabats  from  preceding  Hugoniot 

NOTE:  If  NP[2]=0,  REL  must  immediately  follow  the  call  to  HUG 

for  the  material  of  interest 

1  -  calculate  P  and  release  path  from  input  p  and  E 

2  -  calculate  P  from  input  p,  E.  No  release  path. 

3  -  find  the  lowest  Hugoniot  pressure  that  unloads  to  1  bar, 

i.e.,  the  liquid  side  of  the  steam  dome.  (Subroutine  VAP) 
for  NPF21  =  0 

3  =  0  -  Iterate  on  energy  and  pressure 

1  -  use  E  from  previous  step  like  the  CRALE  codes 

4  =  0  find  closest  point  along  Hugniot  for  release 

1  -  interpolate  along  Hugoniot  to  release  from  exact  values  of  a's 

SETUP 

for  NPT21  =  0 

1  PTOP  Max  P  {Mb}  for  plot  of  Hugoniot  and  adiabats 

2-6  PSAVj  if  PSAV[2]  =  0;  generate  5  adiabats,  incrementing  by  0.2(PTOP) 

or  0.2(ETT);  where  ETT  is  the  value  of  p  at  PTOP 
if  PSAV[2]  *0;  generate  adiabats  starting  at  each  PSAV[2-6]*0. 


forUEI21=1or2 

1  EMUL  value  of  p;  use  p,E  to  calculate  P 

if  material  is  an  HE  and  SETUP(1)=-1  calculate  CJ  point 

2  ENG  value  of  E 

3  constant  o'  added  to  P  for  release  paths 
for  NPf21  =  3 

1  ENG  minimum  value  of  E  for  vapor,  (Es(NE)} 

2  CTE  optional  coefficient  of  thermal  expansion 
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For  NPf21  =  0,  REL  will  generate  release  adiabats  from  selected  pressures  along 
the  Hugoniot  produced  by  the  last  call  to  HUG.  If  SETUP[2]=0.,  SETUP[1]  is  assumed 
to  be  the  peak  stress,  PTOP,  and  unload  paths  are  generated  for  5  stresses  set  to  20, 
40,  60,  80  and  100%  of  PTOP.  If  SETUP[2]*0.,  release  paths  are  generated  starting  at 
stresses  on  the  Hugoniot  corresponding  to  each  value  of  SETUP[2]  through  [6]. 
Repeating  the  line  ROUT=REL  and  NP[2]=0  allows  the  user  to  generate  as  many 
release  paths  as  are  needed. 

For  each  requested  release  stress,  REL  finds  the  data  point  on  the  Hugoniot  that 
equals  or  exceeds  it  and  calls  subroutine  RELSET  to  establish  the  proper  initial 
conditions.  (If  NP[4]  *0,  RELSET  interpolates  in  the  Hugoniot  data  to  start  at  the  exact 
stress  requested.)  Like  HUG,  REL  also  iterates  on  the  energy  density  each  step  to 
converge  on  the  appropriate  point  in  P-V-E  space.  First,  the  top  of  the  cell  is  allowed  to 
move  up  decompressing  the  cell  in  uniaxial  strain.  The  pressure  is  calculated  and  the 
energy  is  then  updated, 

En«En_1+.5(on  +  on.1)dey  (4.6) 

where  the  n-1,n  subscripts  indicate  successive  steps.  The  EOS  is  then  called  again 
with  the  updated  energy,  and  En  is  recomputed  until  it  does  not  change.  As  an  option, 
if  NP[3]  is  set  to  1 ,  REL  skips  the  iteration  and  generates  the  path  exactly  as  would 
CRALE1  and  CRALE2,  i.e.,  using  en  and  En-i.  REL  continues  along  the  release 
'adiabat  untif  either  oy  drops  below  PMN  (1  bar}  or  1 1  is  less  than  EMUMIN  {-.95}. 

\ 

The  printed  output  from  REL  begins  with  a  listing  of  the  data  for  the  point  on  the 
Hugoniot  at  which  it  starts.  REL  generates  a  line-by-line  edit  of  the  unloading  path 
similar  to  that  produced  by  HUG.  At  the  end  of  the  path,  the  initial,  final  and  change  in 
energy  densities  are  printed. 

ROUT=REL  allows  several  special  options  In  addition  to  the  standard  release 
paths  from  the  Hugoniots,  chosen  by  setting  NP[2}>0.  For  NP[2)=  1.  REL  will  produce 
the  unload  adiabat  from  an  arbitrary  input  point  p,  E  set  by  SETUP[1]  and  [2] 
respectively.  REL  calls  EQST  for  P  and  adds  o’  (SETUP[3])  to  get  the  initial  oy  and 
proceeds  to  generate  the  release  adiabat  as  for  NP[2]=0  above. 
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Setting  NPr2V=2.  REL  will  calculate  the  P  for  the  p,  E  input  pair  as  above  but  will 
skip  the  release  path.  This  option  is  used  primarily  to  obtain  single  point  comparisons 
with  other  models  or  data. 

For  NPr21=3.  REL  calls  VAP.  This  routine  starts  by  finding  the  smallest  density 
that  returns  a  1  bar  pressure  for  an  energy  density  equal  to  Em  (the  Pv  cutoff).  This  is, 
presumably,  a  point  on  the  liquid  side  of  the  steam  dome.  It  then  backs  up  the  adiabat, 
by  increasing  p  and  finding  a  P  and  E  that  satisfy  AE=PAV,  until  it  reaches  a  point  that 
also  satisfies  the  Hugoniot  jump  condition,  E  =  -.5P(V-V0),  where  V0  is  the  usual  initial 
specific  volume,  1/p0.  VAP  then  finds  4  additional  release  paths  by  reducing  the  1  bar 
p  by  20%  decrements,  and  increasing  the  E  to  find  another  1  bar  pressure  and  then 
generating  the  new  adiabat.  Currently  these  data  are  only  printed  out,  a  data  set 
suitable  for  plotting  is  not  available. 
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RELEASE  ADIABAT  for  EQS  4  HTH3 
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ENERGY  LOSS  -  init ial , final , ratio  7 . 14 9566E-04  2.039692E-04  2.852889E-01 


4.2.7  PATH  -  Special  Paths. 

The  PATH  subroutine  has  not  been  updated  for  this  version  of  EQS.  Both  the 
routine  and  this  manual  will  be  brought  into  conformity  with  the  rest  of  the  program  in  a 
future  release.  The  following  is  copied  from  a  previous  Users  manual. 

In  addition  to  Hugoniots  and  release  paths,  subroutine  PATH  may  be  used  by 
EQS  to  generate  almost  any  arbitrary  paths  in  stress-strain  space,  including  those  used 
in  the  standard  uniaxial  strain,  hydrostat  and  triax  (i.e.,  uniaxial  stress)  tests.  PATH 
always  produces  a  printout  of  the  path  followed  that  lists  mu,  P,  the  three  principal 
stresses  and  strains,  the  bulk  and  shear  moduli,  yield  surface,  sound  speed,  specific 
energy  and  the  EOS  number.  If  the  energy  is  not  calculated  (the  default  condition)  the 
plastic  strain  is  printed  in  its  place.  The  EOS  number  includes  the  crack  state  as  in 
CRALE1. 

For  uniaxial  load-unload  paths,  two  plots  can  be  produced  by  a  call  to  PATHPL; 
the  first  shows  a  vs  e  and  P  vs  p.  The  second  plots  the  stress  difference  (ax-oy)  vs  P. 
For  triax  runs,  the  stress  difference  is  plotted  vs  the  axial  strain.  A  simple  hydro  run, 
i.e.,  all  three  stresses  equal,  produces  a  plot  with  only  the  P  vs  p  curve. 

The  call  to  PATH  results  in  reading  additional  data  records,  i.e.,  CARD5  . 

,  Usually  the  path  of  interest  requires  several  CARD5s  to  describe  the  full  path. 
^Frequently  parts  of  the  path  are  not  plotted  (e.g.,  the  hydrostatic  loading  which 
precedes  a  triax  run).  Such  segments  of  the  path  may  be  eliminated  from  the  plot  by 
setting  the  appropriate  input  switches. 


1 

l 
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4.2.8  PLOTN-Plotting  Options. 

EQS  saves  the  data  for  each  line  generated  to  allow  for  quite  general  plotting. 
Setting  ROUT=PLOT  will  generate  a  single  frame  containing  one  or  more  data  sets. 

The  axis  of  the  frame  are  determined  by  the  first  line  and  successive  lines  with 
ROUT=PLOT  are  plotted  on  the  same  frame  until  a  line  with  NP[8]=0  is  encountered. 
Setting  NP[8]  to  0  causes  PLOTN  to  advance  the  plot  file  so  that  it  is  ready  to  begin  a 
new  plot.  PLOTN  produces  plots  which  can  overlay  common  K&E  plot  paper.  Hence, 
the  linear  axis  are  multiples  on  1  inch,  the  default  log  scales  correspond  to  the  various 
scales  supplied  by  K&E.  While  this  provides  additional  ease  in  comparing  to  other 
data,  translating  the  plot  files  to  PC  graphics  codes  such  as  DESIGNER  allows  the  user 
to  adjust  the  plot  to  fit  any  format.  The  input  to  PLOT  is; 

m 

1  type  of  plot;  MUST  be  set  to  one  of  the  following: 

2  -  Hugoniot  a  vs  5  from  HUG,  where  8  is  set  by  NP(4) 

3  -  Include  Release  Adiabats  (REL)  with  Hugoniot 

4  -  Material  Parameters,  K,  G,  v,  Cp  and  Us  vs  a  from  HUG 

5  -  a  vs  Up,  from  HUG,  EXP,  or  EQN 

6  -  Us  vs  Up,  from  HUG,  EXP,  or  EQN 
>10  -  Call  PATHPL  to  plot  output  of  PATH 

2  Identifier  of  data  set  to  be  plotted  (Def.  =  last  calc.) 

3-5  see  below 

6  LP:  Line  type,  0=solid  line,  1  <LP<4  for  other  line  types 

7  color  of  data  set 

8  NPP,  continuation  flag:  if  *0,  do  not  advance  plot  frame 

*  this  allows  multiple  data  sets  on  same  plot 

BUT:  NP(8)  on  last  data  set  must  be  0. 

for  NPM1  =  2  or  3  a  vs  p.fi.e.V  or  VA/o 

3  LAX  --  sets  axes  linear  or  logrithmic 

0  -  both  scales  linear 

1  -  Y  (stress)  scale  log 

2  -  both  scales  log 

4  =  LTYP  -  the  parameter  plotted  on  the  X-axis 

0  plot  ay  vs  density,  p 

1  "  "  vs  excess  compression,  p  (i.e.,  AVA/) 

2  "  "  vs  negative  of  volumetric  strain,  -e  (-AVA/o) 

3  "  "  vs  Volume,  V 

,  4  "  "  vs  VA/0 

.  5  If  data  set  is  from  HUG  and 
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NP[5]>0  include  plot  of  P  (E=Eo),  i.e.,  the  hydrostat 
if  data  set  is  from  EXP  and 

NP[5]=0;  plot  different  symbols  for  each  experimenter  and  only  the  name 
of  the  material  in  the  legend 

NP[5]<0;  plot  all  data  with  one  symbol  (=-NP[5])  and  only  the  name 
of  the  material  in  the  legend 

NP[5]=0;  plot  different  symbols  for  each  experimenter  and  list  all 
references  in  the  legend 


for  NP(1)  =  5  or  6 

4  if  >0  plot  release  paths  in  P-Up  space  [NP(1  )=5  only] 

5  if  *0;  plot  references  for  an  EXP  data  set 


SETUP 

for  ALL  plots 

1 

PLXL 

Length  (in  inches)  of  abcissa  (X) 

2 

PLYL 

"  "  "  "  ordinate  (Y) 

3 

PLXMN 

Minimum  value  (in  I/O  units)  of  X  coordinate 

4 

PLXMX 

Maximum .  . . 

5 

PLYMN 

Minimum . .  "Y 

6 

PLYMX 

Maximum . 

In  general,  the  Hugoniot  data  sets  that  are  generated  by  HUG  or  EQN  are  plotted 
as  lines  in  PLOTN,  data  points  from  EXP  are  plotted  by  a  call  to  PLOTH  and  the  output 
of  PATH  are  plotted  by  calling  PATHPL.  The  calls  to  PLOTH  and  PATHPL  are  handled 
by  PLOTN. 


NP£1]=2  produces  a  plot  of  the  principal  Hugoniot,  Figure  4-2,  either  the  line 
generated  by  a  call  to  HUG  or  EQN  or  the  data  set  read  from  the  data  base  by  EXP 
depending  on  the  value  of  NP[2].  The  stress  is  always  plotted  in  the  Y  direction;  the  X 
direction  is  one  of  several  measures  of  compression,  i.e.,  density  (p),  excess 
compression  (p),  total  strain  (e),  volume  (V)  or  specific  volume  (V/V0)  depending  on 
NP[4].  The  axes  may  be  either  linear  or  logrithmic  depending  on  NP[3]  and  up  to  five 
line  types  (NP[6])  and  as  many  as  eight  colors  (NP[7])  are  available  to  help  differentiate 
the  data  sets.  Setting  NPm=3  produces  the  same  plot  as  NP[1]=2  but  includes  the 
release  paths  or  data  as  well. 


When  plotting  the  output  from  HUG  it  is  possible  to  include  a  plot  of  the  zero 
i  energy  hydrostat  by  setting  NP[5]=1 .  Several  options  are  also  available  when  plotting 
'  the  experimental  data  generated  by  EXP.  In  this  case,  if  NP[5]=0,  a  different  symbol  is 
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used  for  each  experimenter's  data  but  only  the  name  of  the  material  is  included  in  the 
legend.  If  NP[5]>0,  each  experimenter  is  also  identified  in  the  legend.  Setting  NP[5]<0 
will  cause  all  of  the  data  for  the  material  to  be  plotted  with  the  same  symbol  as 
determined  by  -NP[5]. 

NP[1]=4  generates  a  plot  of  the  bulk  and  shear  modulus,  poissons  ratio,  and  the 
bulk  and  longitudinal  wave  speeds  calculated  by  HUG  as  functions  of  the  Hugoniot 
mean  pressure,  not  the  stress. 

NP[1]=5  produces  a  plot  of  shock  stress  vs  particle  velocity,  Figure  4-3.  Release 
paths  are  included  if  NP[4]>0  and  the  same  options  on  symbols  and  reference  lists 
based  on  NP[5]  as  for  NP[1]=2  apply.  Line  types  and  colors  are  also  set  as  for  NP[1]=2. 

NP[1]=6  produces  a  plot  of  shock  speed  vs  particle  velocity,  Figure  4-4.  The 
same  options  on  symbols  and  reference  lists  based  on  NP[5]  as  for  NP[1]=2  apply.  Line 
types  and  colors  are  also  set  as  for  NP[1]=2. 


I 
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STRESS  vs  PART  VEL 


0.10  0.20  0.30  0.40  0.50 

Up.  (C/US) 

Figure  4-3.  Example  of  plot  generated  by  ROUT=PLOT  with 
NPfll  =  5. 


SECTION  5 


DATA  BASES 

There  are  two  types  of  data  bases  of  interest  to  the  DNA  user  community,  local 
and  archival.  The  local  data  bases  are  dynamic,  i.e.,  they  reside  on  the  users’  PC  or 
network  and  are  updated  on  a  routine  basis  as  new  information  becomes  available.  The 
current  local  data  bases  include  primary  experimental  shock  and  release  adiabats  and 
material  properties  models.  Archival  data  are  stored  in  DNA  directories  on  the  LANL 
CFS  (Common  File  Storage)  and  generally  are  read-only.  Thus  these  files  are  useful  to 
insure  that  diverse  users  have  access  to  common  data.  Included  in  the  archives  are 
files  relevant  to  both  the  calculational  and  experimental  efforts  conducted  under  the 
HYDROPLUS  yield  verification  program. 

5.1  HUGONIOTS. 

5.1.1  Experimental  Data  (ROCKHUG  and  METHUG). 

The  principal  Hugoniot  is  the  locus  of  points  obtained  by  subjecting  a  material 
initially  at  rest  (ct0,  P0  and  Up  =  0.0)  to  single  shocks  of  differing  magnitudes. 

(Secondary  Hugoniots  start  from  non-zero  initial  conditions.)  Because  a  shock  is  an 
instantaneous  jump  from  one  state  to  another,  there  are  a  set  of  three  equations  which 
relate  the  state  variables,  p[g/cm3],  a[Mb],  E[Te/g],  Us[cm/ps],  and  Up[cm/ps],  i.e. 

a  .  =  p0UsUp  (51) 

pUp  =  [p-pojUs  •  (52) 

E  =  [Up]2/2  <5-3) 

There  are  three  equations  and  five  variables  so  that  any  two  variables  completely 
determine  a  point  on  the  Hugoniot.  Equations  5.1  -  5.3  are  recast  to  relate  the 
dependent  Hugoniot  variables  to  any  pair  in  Table  5-1 .  The  addition  of  one  equation  (or 
system  of  equations,  i.e.  the  EOS  model)  reduces  the  number  of  independent  variables 
to  one.  Thus,  knowledge  of  p  or  p,  the  excess  compression,  coupled  with  an  EOS 
model  is  enough  to  determine  a  unique  point  on  the  Hugoniot  for  comparison  with  data. 
(Note:  As  a  consequence  of  the  jump  conditions,  there  is  an  equal  partioning  of  the 
energy  between  internal  and  kinetic.  As  seen  in  Equation  5.3,  E  and  Up  are  equivalent 
on  the  Hugoniot.) 

^  I 
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Table  5-1 .  Shock  relationships. 


Tergs/g 


U2/  2  [1-Pol J2/G]-1 


[o/p0S]2/2  [1-o/PoS2]-1 


CTp/2p0Tl 


ae/2p0 


[1-PoU2/a]-1 


S/[S-U] 


e=p/ri 


P0U 2 /a 


o/p0S2 


p/tl 


U2/  2 


l/[l-e] 


S/[S-U] 


[o/p0S]2/2  [1-a/PoS2]-1 


[Sp/Til2/2 


[Ss]2/2  1  /[  1  — e] 


U2/2 


a\i  1 2p0r\ 


[Sjj/ti]2/2 


PoU2/a 


U/S 


p/Tl 


Uz/2 


U/S 


o/p0S2 


U2/2 

1/[1-e] 

ae  /  2p0 

[eS]2/2 

if  S=  a  +  b  U  ;  then 


Q 

U 

S 

E 

T\ 

e=p/Tl 

m 

Po[a+bU]U 

— 

a+bU 

U2/2 

a  +  bU 
a  +  U(b-1) 

U 

a  +  bU 

note:  y  -  tn+1]  /  [r|-1]  =  [2-s]  /  e 
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The  principal  Hugoniot  is  very  useful,  since  states  on  it  can  be  obtained 
experimentally  by  measuring  any  two  of  the  state  variables.  A  collection  of  these 
measurements  for  materials  of  interest  to  DNA  resides  on  Units  10  and  1 1 .  These  data 
are  sorted  by  material  in  sets  identified  by  the  author  and  report  ID.  (A  bibliography 
including  all  the  reports  used  in  compiling  the  data  base  is  contained  in  an  ASCII  file 
called  AUTHORS.DAT.)  All  of  the  EOS  parameters  for  a  point  on  the  principal  Hugoniot 
can  be  derived  from  any  pair  of  data  or,  if  the  material's  EOS  is  known,  by  simply 
measuring  a  single  state  variable.  Hence,  each  experimental  Hugoniot  point  is  stored  as 
a  line  of  data  containing  an  identifier  as  to  which  two  variables  are  stored,  the  initial 
density  of  the  sample  and  the  values  of  the  two  variables. 

As  of  1/1/94,  the  materials  included  in  the  Unit  10  database  (rocks,  etc)  are  listed 
in  Table  5-2;  those  on  Unit  1 1  (metals,  plastics,  etc)  in  Table  5-3.  Table  5-4  contains  a 
partial  listing  of  the  data  for  dry  calcite  and  is  an  example  of  the  data  stored  on  these 
databases 

The  DNA  shock  database  includes  data  stored  in  LLL  and  LANL  data  libraries, 
either  referenced  to  the  original  author  or  by  inference  to  the  library  from  which  it  was 
obtained.  While  our  primary  source  is  the  original  author's  report,  when  these  are  not 
available  the  data  in  the  DNA  file  is  simply  a  copy  of  the  LLL  or  LANL  library 

5.1.2  Shock  Velocity,  Us,  vs  Particle  Velocity, Up,  Fits  (USUP.dat,  1-1-94). 

It  has  been  observed  that  for  many  materials,  over  a  wide  range  of  velocities,  the 
%  shock  velocity  increases  linearly  with  particle  velocity,  i.e.; 

Us  =  C0  +  sUp  (5-4) 

where  C0  and  s  are  found  by  least-squares  fits  to  the  data.  With  these  two  coefficients 
any  point  on  the  Hugoniot  is  determined  by  specifying  only  one  variable.  Because  of 
this  simplification,  many  experimenters  and  modelers  use  Equation  5.4  to  fit  data  or 
develop  material  models.  Some  of  these  fits  have  been  collected  in  a  file,  Table  5-5 
and  USUP.dat,  (Unit  14),  which  can  be  accessed  via  EQN  in  the  driver  program 
(Section  4)  to  generate  a  complete  Hugoniot  for  plotting  and  comparisons  with  other 
lines.  Please  note  that  Equation  5.4  should  not  be  used  when  the  pressure  range 
(  includes  phase  transitions. 
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Table  5-2.  Materials  stored  in  the  ROCKHUG.dat  file  (Unit  10)  11-11  -93. 


in  order  on  file 

in  numerical  order 

WATER 

108. 

WATER 

108. 

ICE 

108.1 

ICE 

108.1 

AL2O3/EPOXY(40%) 

400.1 

DRY  ALLUVIUM 

210.1 

DRY  ALLUVIUM 

210.1 

WET  ALLUVIUM 

210.2 

WET  ALLUVIUM 

210.2 

QUARTZ 

220. 

ANORTHOSITE 

221. 

WET  QUARTZ 

220.2 

BASALT 

260. 

POROUS  QUARTZ 

220.3 

ME  GROUT 

330. 

FUSED  QUARTZ 

220.4 

GROUT 

330.1 

ANORTHOSITE 

221. 

CONCRETE 

330.3 

GRANITE 

222. 

RHYOLITE 

330.2 

SANDSTONE 

223. 

GROUT 

330.3 

TONALITE 

224. 

GNEISS 

370. 

DRY  TUFF 

230.1 

DRY  CALCITES 

240. 

WET  TUFF 

230.2 

SINGLE  CRYSTAL  CALCITE 

240.1 

DRY  CALCITES 

240. 

WET  CALCITES 

240.2 

SINGLE  CRYSTAL  CALCITE 

240.1 

ARAGONITE  (hi-P  CALCITE) 

242. 

WET  CALCITES 

240.2 

MISC  ROCK  (RHO<3.0) 

310.1 

ARAGONITE  (hi-P  CALCITE) 

242. 

MISC  ROCK  (RHO>3.0) 

310.2 

BASALT 

260. 

MISC  ROCK  (RHO>4.0) 

310.3 

SALT 

270. 

PUMICE 

331. 

SHALE 

280. 

GRANITE 

222. 

PHYLLITE 

280.1 

QUARTZ 

220. 

MISC  ROCK  (RHO<3.0) 

310.1 

WET  QUARTZ 

220.2 

MISC  ROCK  (RHO>3.0) 

310.2 

POROUS  QUARTZ 

220.3 

MISC  ROCK  (RHO>4.0) 

310.3 

FUSED  QUARTZ 

220.4 

ME  GROUT 

330. 

SANDSTONE 

223. 

GROUT 

330.1 

SALT 

270. 

CONCRETE 

330.3 

SHALE 

280. 

RHYOLITE 

330.2 

PHYLLITE 

280.1 

GROUT 

330.3 

TONALITE 

224. 

PUMICE 

331. 

DRY  TUFF 

230.1 

GNEISS 

370. 

WET  TUFF 

230.2 

AL2O3/EPOXY(40%) 

400.1 
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Table  5-3.  Materials  stored  in  the  METHUG.dat  file  (Unit  11)  11-11  -93. 

in  order  on  file  in  numerical  order 


ALUMINUM -AL  13. 

ANTIMONY -SB  51. 

BERYLLIUM  -  BE  4. 

BISMUTH  -  Bl  83. 

CADMIUM  -  CD  48. 

CHROMIUM  -  CR  24. 

COBALT  -  CO  27. 

COPPER  -  CU  29. 

GOLD  -  AU  79. 

INDIUM  -  IN  49. 

IRON  -  FE  (LLL  NO.  41)  26. 

LEAD  -  PB  82. 

MAGNESIUM  -  MG  12. 

MOLYBDENUM-  MO  42. 

NICKEL-  Nl  28. 

NIOBIUM -NB  41. 

PALLADIUM  PD  46. 

PLATINUM  -  PT  78. 

RHODIUM  -  RH  45. 

SILVER  -  AU  47. 

TANTALUM  -  TA  73. 

THALLIUM -TL  81. 

THORIUM  -  TH  90. 

TIN  -  SN  50. 

TITANIUM  -  Tl  22. 

TUNGSTEN  -  W  74. 

WC  -  TUNGSTEN  CARBIDE  74.1 
VALLADIUM  -  V  23. 

ZINC  -  ZN  30. 

ZIRCONIUM  -  ZR  40. 

BRASS  29.5 

LITHIUM  DEUTERIDE  3.2 

LITHIUM  FLORIDE  3.1 

PMMA  12.5 

TEFLON  12.6 


LITHIUM  FLORIDE 

3.1 

LITHIUM  DEUTERIDE 

3.2 

BERYLLIUM  -  BE 

4. 

MAGNESIUM  -  MG 

12. 

PMMA 

12.5 

TEFLON 

12.6 

ALUMINUM  -  AL 

13. 

TITANIUM  -  Tl 

22. 

VALLADIUM  -  V 

23. 

CHROMIUM  -  CR 

24. 

IRON-FE  (LLL  NO.  41) 

26. 

COBALT  -  CO 

27. 

NICKEL-  Nl 

28. 

COPPER  - CU 

29. 

BRASS 

29.5 

ZINC  -  ZN 

30. 

ZIRCONIUM  -  ZR 

40. 

NIOBIUM  -  NB 

41. 

MOLYBDENUM-  MO 

42. 

RHODIUM  -  RH 

45. 

PALLADIUM  PD 

46. 

SILVER  -  AU 

47. 

CADMIUM  -  CD 

48. 

INDIUM  -  IN 

49. 

TIN  -  SN 

50. 

ANTIMONY -SB 

51. 

TANTALUM  -  TA 

73. 

TUNGSTEN  - W 

74. 

WC  -  TUNGSTEN  CARBIDE 

74.1 

PLATINUM  -  PT 

78. 

GOLD  -  AU 

79. 

THALLIUM  -  TL 

81. 

LEAD  -  PB 

82. 

BISMUTH  -  Bl 

83. 

THORIUM  -  TH 

90. 

I 

I 
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Table  5-4.  Example  of  experimental  Hugoniot  data  stored 
on  the  ROCKHUG.dat  and  METHUG.dat  files. 


DRY  CALCITES 

HE  16.  5.90 

ID  2.29 

SP  3.02  5.50 

3554 

KTECH  (92e)  SALEM  LIME  9-18 

SP  2.40 

8.80 

3556 

SP  2.35 

20.6 

3558 

SP  2.52 

35.7 

3561 

SP  2.26 

2.00 

10.3 

3564 

ID  2.7 

SU  5.15 

.86 

.5' 

COLEMAN(92)  DM1 

SU  4.83 

.805 

1' 

SU  4.83 

.74 

1' 

SU  4.55 

.6 

2' 

SU  4.47 

.595 

2.5' 

ID  2.70 

SP  6.094 

11.84 

3467 

KTECH  (92a)  DM  MARBLE  GAFFNEY  4- 

SP  4.797 

14.3 

3469 

PS  4.0 

54.3 

M 

UP  0.91 

102.5 

3527 

UP  2.69 

1.12 

156.4 

3513 

ID  2.70 

UP  .051 

12.4 

3470 

KTECH  (92a)  DM  MARBLE  GAFFNEY  4. 

UP  .571 

48. 

3471 

UP  .569 

50. 

3473 

UP  .054 

12.5 

3474 

UP  .532 

51. 

3475 

UP  .555 

50. 

3476 

ID 

UP  2.8  .118 

20.1 

3489 

KTECH  (92f)  LiME 

UP  2.79 

.119 

18.3 

3490 

UP  2.71 

.200 

24.3 

3499 

UP  2.69 

.223 

26.1 

3492 

UP  2.68 

.220 

28.7 

3491 

UP  2.75 

.372 

54.2 

3497 

UP  2.74 

.377 

59.7 

3498 

UP  2.73 

.293 

34.4 

3503 

UP  2.73 

.345 

39.2 

3505 

UP  2.75 

.415 

54.7 

3504 

ID 

UP  2.694  0.11415.4 

DM6 

FURNISH(92A)  DANBY  MARBLE  (4-17) 

RA 

UP  2.694  0.605  66.5 

DM7 

DM6 

RA 

UP  2.694  0.679  83.7 

DM8 

DM7 

ND 

****  end  of  DRY  CALCITES  data 

64 


Table  5-5.  Material  Us-Up  relations  available  on  USUP.dat  file. 


ID  Num 

Density 

(g/cc) 

Material 

Authors 

108. 

1.0 

Water 

CRT  (92) 

108.31 

1.003 

Water 

BAKANOVA(76) 

108.1 

0.91 

Ice 

CRT  (92) 

211.31 

2.21 

Clay 

TRUNIN(88) 

222.1 

2.62 

Granite 

Aherns  et  al 

222.31 

2.60 

Granite 

TRUNIN(88) 

223.31 

Sandstone 

TRUNIN(88) 

224.31 

2.76 

Siltstone 

TRUNIN(88) 

231.31 

2.74 

Tuff 

TRUNIN(88) 

231.32 

1.89 

DZ  Tuff 

CRT  (92) 

240.31 

2.71 

Danby  Marble 

CRT  (92) 

240.31 

2.72 

Limestone 

TRUNIN(88) 

240.32 

2.84 

Dolomite 

TRUNIN(88) 

310.21 

2.74 

Syenites 

TRUNIN(88) 

310.32 

2.89 

Gabbro 

TRUNIN(88) 

310.33 

2.62 

Porphyrites 

TRUNIN(88) 

310.34 

Igneous  rock 

TRUNIN(88) 

310.35 

330.0 

1.95 

it 

MISTY  ECHO  Grout 

CRT  (92) 

330.10 

1.99 

MISTY  ECHO  Grout 

GERMAIN(92) 

330.11 

1.937 

HPNS12  Grout 

GERMAIN(92) 

333.31 

2.77 

Shales 

TRUNIN(88) 

1. 

1.193 

LEXAN 

Aherns  et  al 

1. 

0.055 

Polystyrene 

Aherns  et  al 

3.130 

2.641 

LiF 

FURNISH  Aug(90) 

3.140 

2.64  . 

LiF 

STEINBERG  (91),  AHERNS, 

13.0 

2.785 

Aluminum  2024 

AHERNS  ET  AL  (OCT  91) 

13.31 

Aluminum 

TRUNIN(88) 

13.32 

2.78 

Aluminum  2024L 

LANL(80) 

13.33 

2.70 

Aluminum  6061 L 

LANL(80) 

13.40 

2.70 

Aluminum  6061 S 

STEINBERG(91) 

13.35 

2.80 

Aluminum  7075L 

LANL(80) 

13.30 

2.689 

Aluminum 

FURNISH  Aug(90) 

29.31 

Copper 

TRUNIN(88) 

26.0 

7.86 

Steel  4340 

CRT(92) 

74.1 

14.9 

WC 

CRT  fit  for  KTech  Tests 

74.11 

14.9 

WC 

CRT(92) 

74.1 

14.9 

WC 

STEINBERG  (91) 

73.30 

16.66 

Tantalum 

AHERNS  ET  AL 

LANL(80) 
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5.2  RELEASE  ADIABATS  (RELAD.dat,  1-1-94). 


The  data  describing  the  release  paths  are  not  stored  as  conveniently  as  the 
Hugoniot  input  since  unknown  values  cannot  be  calculated  from  the  measured  data. 
Rather  the  suite  of  four  variables  (a,  p,  Up,  and  E)  at  points  along  a  release  path  must 
be  derived  from  the  measurements  and  stored  in  the  data  base.  Both  1  d  code 
calculations  and  LASS  analysis  are  used  to  obtain  the  release  paths. 

Recent  experiments  in  which  1 D  codes  are  used  to  derive  the  complete  set  of 
material  properties  are  currently  being  incorporated  into  EXP  and  will  become  available 
in  a  later  release.  The  current  [1/1/94]  release  adiabat  database  ( RELAD.dat  on  Unit 
12)  contains  the  105  paths  listed  in  Table  5-7.  All  of  these  data  were  obtaqined  in 
reverse  ballistic  tests  conducted  by  SNL;  the  shot  names  correspond  to  the  Hugoniot 
data  stored  in  the  ROCKHUG  and  METHUG  data  files. 


An  example  of  the  data  for  a  single  release  path  is  included  in  Table  5-6.  The 
data  are  stored  in  ASCII  format  to  allow  them  to  be  easily  updated  and  ported  to  any 
PC  or  workstation. 


Table  5-6.  ExampleofdataonRELAD.dat  (1-1-94). 


\ 


I 


P(MB) 

23  BUL6 

4.1033E-01 

4.0463E-01 

4.031 9E-01 

4.0221  E-01 

4.0196E-01 

4.0029E-01 

3.9900E-01 

3.9803E-01 

3.9755E-01 

3.9100E-01 

3.6757E-01 

3.461  IE-01 

3.2386E-01 

3.0734E-01 

2.91 51  E-01 

2.7687E-01 

2.6395E-01 

2.5293E-01 

2.4484E-01 

2.41 81  E-01 

2.3527E-01 

2.2830E-01 


DEN(g/cc) 

22 

3.8290E+00 
3.8159E+00 
3.8142E+00 
3.8135E+00 
3.8133E+00 
3.8109E+00 
3.8087E+00 
3.8077E+00 
3.8069E+00 
3.7983E+00 
3.7631  E+00 
3.7269E+00 
3.6847E+00 
3.6497E+00 
3.6128E+00 
3.5753E+00 
3.5391  E+00 
3.5055E+00 
3.4791  E+00 
3.4688E+00 
3.4458E+00 
3.4199E+00 


Up(cm/us) 

3.361 0E-01 

3.3837E-01 

3.3877E-01 

3.3898E-01 

3.3905E-01 

3.3957E-01 

3.4001  E-01 

3.4027E-01 

3.4044E-01 

3.4241  E-01 

3.5000E-01 

3.5745E-01 

3.6571  E-01 

3.7227E-01 

3.7893E-01 

3.8545E-01 

3.9153E-01 

3.9699E-01 

4.01 17E-01 

4.0278E-01 

4.0633E-01  ' 

4.1025E-01 


E(Terg/g) 

5.6485E-02 
5.6124E-02 
5.6078E-02 
5.6060E-02 
5.6052E-02 
5.5986E-02 
5.5926E-02 
5.5899E-02 
5.5876E-02 
5.5648E-02 
5.4770E-02 
5.3905E-02 
5.2944E-02 
5.2166E-02 
5.1372E-02 
5.0589E-02 
4.9852E-02 
4.9183E-02 
4.8662E-02 
4.8457E-02 
4.8009E-02 
4.751 5E-02 
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Table  5-7.  List  of  105  sets  of  release  Adiabats  on  RELAD.dat  (1-1-94). 


SHOT  NAME 

num 

Phug 

SHOT  NAME 

num 

Phug 

(mB) 

(mB) 

1 

DNA1 

50 

0.0419 

51 

HT4 

50 

0.0419 

2 

DNA2 

50 

0.0617 

52 

HT5 

50 

0.0851 

3 

DNA3 

50 

0.0834 

53 

HT6 

50 

0.0161 

4 

DNA4 

50 

0.1024 

54 

HT7 

50 

0.1068 

5 

DNA5 

50 

0.1197 

55 

HT8 

50 

0.1397 

6 

DNA6 

50 

0.1187 

56 

HT9 

32 

0.3504 

7 

DNA7 

50 

0.1512 

57 

HT10 

20 

0.5343 

8 

DNA8 

50 

0.1854 

58 

HT11 

25 

0.5717 

9 

DNA9 

19 

0.0487 

59 

HT12 

50 

0.0357 

10 

DNA11 

20 

0.0211 

60 

HT13 

50 

0.0770 

11 

DNA12 

12 

0.1639 

61 

BEX1 

35 

0.0996 

12 

DNA13 

22 

0.2457 

62 

BEX2 

50 

0.0385 

13 

DNA14 

20 

0.2861 

63 

BEX3 

31 

0.0646 

14 

DNA15A 

47 

0.3519 

64 

BEX4 

36 

0.1449 

15 

DNA16 

19 

0.4356 

65 

BEX5 

49 

0.1160 

16 

DNA18 

15 

0.4554 

66 

BEX6 

50 

0.0502 

17 

DNA19 

14 

0.5542 

67 

BEX7 

50 

0.0923 

18 

DNA20 

15 

0.6592 

68 

BEX8 

50 

0.1568 

19 

BUL2 

50 

0.0225 

69 

BEX11 

50 

0.1222 

20 

BUL3 

43 

0.1122 

70 

BEX13 

25 

0.1717 

21 

BUL4 

28 

0.2751 

71 

DLS1 

36 

0.0109 

22 

BUL5 

16 

0.4689 

72 

DLS3 

30 

0.0320 

23 

BUL6 

22 

0.4103 

73 

DLS6 

26 

0.1218 

24 

BUL7 

16 

1.0117 

74 

DLS8 

34 

0.0576 

25 

BUL8 

11 

0.1750 

75 

ILS1 

50 

0.0570 

26 

BUL 

16 

0.0503 

76 

ILS2 

41 

0.0792 

27 

NT2 

42 

0.0188 

77 

ILS3 

34 

0.0954 

28 

NT3 

50 

0.0640 

78 

ILS9 

50 

0.0129 

29 

NT4 

50 

0.0973 

79 

ILS10 

50 

0.0252 

30 

NT5 

50 

0.1196 

80 

1LS12 

50 

0.0449 

31 

NT6 

48 

0.0428 

81 

ILS13 

47 

0.0060 

32 

PT3 

50 

0.0577 

82 

ILS14 

48 

0.0405 

33 

PT4 

44 

0.0873 

83 

ILS15 

45 

0.0943 

34 

PT5 

50 

0.1329 

84 

PF2 

50 

0.0421 

35 

•PT6 

46 

0.0181 

85 

DM1 

49 

0.0325 

36 

RT1 

50 

0.0155 

86 

DM2 

50 

0.0367 

37 

RT3 

*  27 

0.0603 

87 

DM3 

50 

0.0561 

38 

RT4A 

50 

0.1196 

88 

DM4 

50 

0.0985 

39 

RT5 

50 

0.0867 

89 

DM5 

40 

0.1302 

40 

RT6 

50 

0.0220 

90 

DM6 

50 

0.0160 

41 

TAN2 

50 

0.0778 

91 

DM7 

50 

0.0664 

42 

TAN3 

50 

0.1364 

92 

DM8 

50 

0.0830 

43 

TAN4 

50 

0.2299 

93 

DM9 

50 

0.1303 

44 

TAN7 

17 

0.4485 

94 

DM14 

50 

0.0443 

45 

TAN9 

50 

0.0293 

95 

DM15 

50 

0.0463 

46 

DZ1 

50 

0.0476 

96 

DM16 

50 

0.0586 

47 

DZ2 

25 

0.1249 

97 

DM17 

50 

0.0847 

48 

DZ3 

10 

0.3148 

98 

HT16 

31 

0.25817 

49 

DZ4 

13 

0.3866 

99 

HT17 

31 

0.20502 

50 

DZ5 

21 

0.5494 

100 

SLP1 

40 

0.05232 

101 

SLP2 

40 

0.08612 

102 

SLP3 

30 

0.16274 

103 

SLP4 

27 

0.26580 

104 

SLP6 

11 

1.39710 

105 

SLP7 

26 

0.65713 
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5.3  MATERIAL  MODELS  (EOSS.dat)  1-1-94. 


In  addition  to  the  experimental  data  bases,  TRT  maintains  a  file  containing  the 
input  decks  for  a  large  number  of  material  models,  for  example  there  are  24  different  HE 
materials  available.  While  the  emphasis  of  the  table  was  to  preserve  the  models  of 
various  geologic  materials  used  in  DNA  calculational  efforts,  EOS  models  formetals, 
plastics  and  other  materials  that  were  developed  under  DNA  and  WES  programs  are 
included  in  the  data  file.  The  materials  whose  EOS  models  are  store  in  the  EOSS.dat 
data  base  are  listed  in  Table  5-8  below. 


Table  5-8.  Listing  of  the  EOS  models  in  the  EOSS.dat  data  file, 
Tape  14  as  of  10/27/93. 


rocks  -  water  --  ice  -  grouts 
ICE  H20N  HTH3  DZ2 

LMS023  LMST01  DM2GR1  DMM001 

SHALE  BXR6  BXR8  BXR12 


DZ5  TEN6  TUF58 

LTRM20  SBAGS 

BXR13  MES2  MES195 


metals  ~  plastics  -  etc 
STEEL  AL  WC  TA 


TAHY  TAHHY  TEFLON  FOAM40  PMMA  LIF 


High  Explosives 

AIR 

ANFO 

ATX 

H6 

MEN2 

NM 

PETN 

PETNA 

PETNB 

RX39A 

RX39B 

TNT85 

COMPB 

NM78 

PETNC 

TNT83 


C4 

PBX944 

PETNS 

TNT72 


DSHTA 

PXN103 

9404 


DSHTC 

N103 
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5.4  HYDROPLUS  ARCHIVAL  DATABASES. 

5.4.1  Introduction. 

To  assist  in  maintaining  the  documentation  required  by  the  HYDROPLUS  yield 
verification  program,  a  single  collection  point  for  data  was  established  for  DNA  on  the 
CFS  (Common  File  Storage)  system  at  LANL.  A  number  of  different  organizations  and 
agencies  were  involved  in  fielding  gages,  measuring  equation  of  state  properties, 
calculating  probable  ground  response  and  making  yield  estimates.  The  nature  of  the 
HYDROPLUS  yield  estimation  required  close  communication  between  these  groups  and 
the  ability  to  store  and  share  files.  The  Common  File  System  (CFS)  at  LANL  was 
chosen  as  a  suitable  location  for  such  a  database  since  all  of  the  organizations  have 
access  to  it.  Titan  Research  and  Technology,  TRT  (formerly  CRT),  created  two 
directories  to  serve  as  the  repositories  on  the  CFS.  The  first  is  primarily  for  equation  of 
state  data  and  numerical  calculations  and  is  unclassified;  the  second,  located  in  the 
classified  partition  of  the  CFS,  contains  summaries  of  the  data  from  the  five 
HYDROPLUS  nuclear  tests  fired  to  date. 

5.4.2  Equation  of  State/Numerical  Results. 

The  EOS  database  contains  Hugoniot  points  and  release  adiabat  paths  for 
specific  materials  of  interest  to  the  HYDROPLUS  program.  These  have  been  culled 
from  the  open  literature,  i.e.,  the  files  described  in  Section  5.1,  and  augmented  by 
specific  HYDROPLUS-sponsored  high-velocity  impact  tests  on  small  samples.  The 
materials  included  are  tabulated  below.  Release  adiabats  for  materials  tested  by 
x  Sandia  are  contained'  in  a  single  file,  rel-sum.std.  These  materials  include  tuff,  rhyolite, 
and  limestone.  A  complete  list  of  the  files  in  the  unclassified  database  is  included  in 
Appendix  A  Several  text  files  [readme.--]  are  interspersed  throughout  the  file  to  aid  the 

user  in  finding  any  specific  data  set. 

In  addition  to  the  experimental  EOS  data  sets  listed  in  Table  5-9,  calculational 
models  are  also  included  in  the  database.  Each  of  the  three  organizations  that 
participated  in  HYDROPLUS  ground  shock  calculations  [CRT,  SAIC.  and  SCUBED]  has 
a  subdirectory  in  which  to  collect  their  relevant  files.  All  users  have  access  rights  which 
allow  them  to  add  files  to  the  notredame  directory,  however,  once  the  files  have  been 
created,  only  TRT  can  delete  or  modify  them.  This  scheme  is  designed  to  centralize 
1  changes  and  maintain  accountability. 
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Table  5-9.  Materials  in  HYDROPLUS  EOS  database. 


Material 

file(s) 

date  of  last  entry 

rhyolite 

bexarhug.std,  rhyolite.std 

1-22-91,  3-8-93 

wet  tuff 

bullionhug.std, 
tenabohug.std,  wet  tuff.std 

1-25-91,  1-25-91,  3-8-93 

dry  tuff 

dry_tuff.std 

3-8-93 

ME  grout 

mistyechohug.std, 

megrout.std 

1-25-91,  3-8-93 

low-porosity  limestone, 
marble,  calcite,  aragonite 

sol-lime-plus.s 

3-8-93 

dry  calcite 

dry^calcite.std 

3-8-93 

wet  calcite 

wet__calcite.std 

3-8-93 

single  crystal  calcite 

sngcryscalcite.s 

3-8-93 

shale 

shale.std 

3-8-93 

water 

water,  std 

3-8-93 

ice 

ice.std 

3-8-93 

AI203/epoxy(40%) 

epoxy,  std 

3-8-93 

misc.  grout 

grout,  std 

3-8-93 

sandstone 

sandstone.std 

3-8-93 

These  files  are  all  stored  at  Los  Alamos  (LANL)  on  the  common  file  system 
(CFS)  in  the  root  directory  /notredame.  Potential  users  may  contact  DNAFC,  Eric 
Rinehart,  for  the  password.  The  files  are  predominantly  stored  in  standard  format. 
Standard  format  is  produced  by  running  STEXT  on  an  ASCII  file.  The  file  may  then  be 
,  retrieved  onto  a  variety  of  machines  using  a  local  NTEXT  utility  which  converts  the  files 
vback  into  native  ASCII  format  on  the  users'  machine.  This  makes  is  possible  to  retrieve 
the  files  on  such  diverse  systems  as  a  VAX  or  UNIX  machine.  Once  there,  any 
standard  text  editor  or  word  processor  may  be  used  to  view  the  file. 


I 
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5.4.3  Nuclear  Shot  Data. 

There  have  been  five  nuclear  tests  on  which  HYDROPLUS  gages  were  fielded 
and  a  yield  estimate  made.  BULLION  and  TENABO  were  'standard'  tests  as  regards 
the  Threshhold  Test  BanTreaty  (TTBT).  Both  events  were  shots  fired  in  canisters 
emplaced  in  vertical  drill  holes  in  tuff  at  NTS.  BEXAR  was  emplaced  in  a  vertical  drill 
hole  in  rhyolite.  DISTANT  ZENITH  and  HUNTERS  TROPHY  were  horizontal  tunnel 

shots  in  tuff. 

The  purpose  of  maintaining  a  HYDROPLUS  database  is  to  provide  the 
HYDROPLUS  community  with  a  single,  definitive  set  of  results  for  each  event  in  which  it 
participates.  The  bibliography  accompanying  the  data  provides  the  necessary  paper 
trail  and  a  library  of  the  referenced  reports  is  maintained  at  TRT. 

Each  nuclear  test  has  been  summarized  in  an  EXCEL  spreadsheet.  This  allows 
the  inclusion  of  graphics  to  illustrate  the  configuration  of  each  test.  While  each 
summary  has  a  slightly  different  format  to  accomodate  the  unique  configuration  of  each 
test,  in  general  there  are  sections  describing: 

•  geometry 

•  media 

•  yield 

•  gage  locations 

.  .gage  measurements,  both  the  initial  quick  look  and  latest  best  value  and  the 
free-field  equivalent,  i.e.,  borehole  corrected 

•  plots  of  peak  pressure  and  velocity  vs  range 

•  references 

These  spreadsheets  are  stored  on  CFS  in  the  directory  /hat.  The  formats  are. 

1 .  The  encoded  binary  EXCEL  file  (.EXL).  This  file  is  encoded  with  UUENCODE;  a 
copy  of  which  is  available  to  users  in  the  unclassified  /notredame  directory,  Appendix 
A.  These  are  shareware  routines  that  preserve  the  PC  file  format. 

2.  An  ASCII  equivalent  text  file  (.TXT)  suitable  for  viewing/editing  on  the  DNA 
CRAY/VAX  system  at  LANL.  The  graphics  cannot  be  included  in  this  file. 

3.  A  comma-delimited  version  (.CSV)  for  reading  into  spreadsheet  programs  that  can 
not  rcaJ  EXCEL  files.  As  might  be  expected,  This, format  does  not  contain  graphics. 

'  An  unclassified  example  of  the  HUNTERS  TROPHY  data  is  included  as  Appendix  B. 
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LIST  OF  FILES  IN  UNCLASSIFIED  DIRECTORY 
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EXAMPLE  OF  THE  HYDROPLUS  NUCLEAR  SHOT  DATA  BASE 
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